PRECIPITACIÓN EN PUEBLA
Precipitación v4:
https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network-monthly
https://www.ncei.noaa.gov/data/ghcnm/v4/precipitation/
https://www.ncei.noaa.gov/pub/data/ghcn/v4/products/StationPlots/MX/
Nombres:
https://www.ncei.noaa.gov/pub/data/ghcn/v4/
https://www.ncei.noaa.gov/data/ghcnm/v4/precipitation/doc/
Ubicación de la estación: 19.0000, -98.1833
CONFIGURACIÓN DE LA NOTETBOOK¶
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="whitegrid")
# The color palette is made up of the 20 colors. Hex color codes: #395e77, #413344, #614c65, #806485, #936397, #a662a8, #664972, #463c57, #6e8da9, #91bcdd, #567d99, #305662, #264d4d, #315c45, #8a9a65, #b6b975, #b65d54, #b60033, #98062d and #800022.
# https://colorkit.co/palette/413344-614c65-806485-936397-a662a8-664972-463c57-6e8da9-91bcdd-567d99-395e77-305662-264d4d-315c45-8a9a65-b6b975-b65d54-b60033-98062d-800022/
sns.set_palette(["#194a7a", "#b60033", "#315c45", "#b65d54", "#B9445F", "#567d99", "#395e77", "#413344", "#614c65", "#806485", "#936397", "#a662a8", "#664972", "#463c57", "#6e8da9", "#91bcdd", "#305662", "#264d4d", "#8a9a65", "#b6b975", "#98062d", "#800022"])
sns.color_palette()
# Personalización global con matplotlib
plt.rcParams.update({
'axes.titlesize': 14, # Tamaño del título
'axes.titleweight': 'bold', # Negrita en el título
'xtick.labelsize': 8, # Tamaño de los xticks
'ytick.labelsize': 8, # Tamaño de los yticks
'grid.color': 'gray', # Color de las líneas del grid
'grid.linestyle': '--', # Estilo de línea (puede ser '-', '--', '-.', ':')
'grid.linewidth': 0.5, # Grosor del grid
'axes.grid': True, # Asegura que el grid esté activado
'axes.grid.axis': 'both', # Aplica el grid a ambos ejes
'lines.linewidth': 1.2, # Grosor de las líneas
'figure.figsize': (12, 6), # Tamaño de la figura
})
CARGAR DATOS¶
data=pd.read_csv('MXN00021035.csv')
pre=data.iloc[:,6] # Precipitacion, es la columna 5
date=data.iloc[:,5] # Date, es la columna 6
date = date.astype(str).str.replace(r'(\d{4})(\d{2})', r'\1/\2', regex=True) # La fecha está como 195210 y la pasamos a 1952/10
date = pd.to_datetime(date, format='%Y/%m') # Lo convertimos en fecha
pre = pd.Series(pre.values, index=date) # Creamos una Serie
pre
195209
1952-10-01 174
1952-11-01 592
1952-12-01 0
1953-04-01 47
1953-05-01 137
...
2009-08-01 1474
2009-09-01 3161
2009-10-01 1190
2009-11-01 155
2009-12-01 51
Length: 672, dtype: int64
import matplotlib.dates as mdates
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3)) # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) # Formato de fecha
plt.plot(pre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla")
plt.savefig('imagenes/01-01-precipitacion.svg', bbox_inches='tight')
plt.show()
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12)) # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) # Formato de fecha
plt.plot(pre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla (2000-2010)")
plt.xlim(pd.Timestamp('1990-01-01'), pd.Timestamp('2010-01-01'))
plt.savefig('imagenes/01-02-precipitacion-zoom.svg', bbox_inches='tight')
plt.show()
TRAIN Y TEST¶
# Partir la serie para train y test
pre_total = pre.copy() # Copia de la serie original
# Todas hasta los ultimos 12 meses
pre = pre_total[:-12] # Entrenamiento: todos menos los últimos 12 meses
pre_test = pre_total[-12:] # Test: últimos 12 meses
SERIE SIN TRANSFORMACIONES, SOLO ESTANDIARIZADA¶
Como la vamos a comparar con la transformacion de Yeo-Johnson, estandarizamos para comparar los modelos.
mu = pre.mean() # Media
sigma = pre.std() # Desviación estándar
# Normalizar la serie
zpre = (pre - mu) / sigma # Estandarizar la serie
ESTACIONARIEDAD¶
from statsmodels.tsa.stattools import adfuller
adfuller(zpre)
(-6.150368236022414,
7.588129761677003e-08,
11,
648,
{'1%': -3.4404817800778034,
'5%': -2.866010569916275,
'10%': -2.569150763698369},
1333.5800521506967)
Estadístico ADF = -6.150
Valor-p = 7.59e-08
Número de rezagos = 11
Número de observaciones = 648
Valores críticos:
1% -> -3.4405
5% -> -2.8660
10% -> -2.5692
Log-likelihood = 9999.935
- El estadístico ADF es $-6.150$, menor que los valores críticos a los niveles del 1%.
- El valor-p es $7.59 × 10^{-8}$ mucho menor al nivel de significancia $\alpha$ = 0.05.
Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.
AUTOCORRELACIONES¶
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import * # Este es el archivo que se usa para el curso
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(zpre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/02-01-fac-facp-serie.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(zpre), acf(zpre, nlags=48)[1:] )
Valores de autocorrelacion significativos: r1: 0.5114347926277153 r2: 0.23788150833700863 r4: -0.31136216826642676 r5: -0.5289967604156934 r6: -0.5765988530134252 r7: -0.4540298620448335 r8: -0.24420551667402526 r10: 0.2960004737880804 r11: 0.5240658344251691 r12: 0.6173499572874946 r13: 0.424401244621075 r14: 0.18265000690415573 r16: -0.31903920298321264 r17: -0.514286815393446 r18: -0.527829591696689 r19: -0.4187324328138189 r22: 0.30067304144252044 r23: 0.5028304240974979 r24: 0.5626631891391944 r25: 0.3994834359095428 r28: -0.34516091598549786 r29: -0.48547226384609826 r30: -0.5033566062877032 r31: -0.38761684940293784 r34: 0.3085645527226526 r35: 0.4965812730150526 r36: 0.553098035132298 r37: 0.36241966295184613 r40: -0.3454697704347722 r41: -0.4674618422218585 r42: -0.46761339396960605 r43: -0.34243418620086286 r46: 0.3414843279972229 r47: 0.48860016886221763 r48: 0.5000880408130988
facp = FACP(len(zpre), pacf(zpre, nlags=48)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.5122108697030228 rho 3: -0.1371024981440388 rho 4: -0.35568794886746186 rho 5: -0.3390668723583037 rho 6: -0.26990992305506756 rho 7: -0.13651941979413404 rho 8: -0.082068919602615 rho 11: 0.19446853660047775 rho 12: 0.25021094781029624 rho 19: -0.0924235210146644 rho 24: 0.12444505947195694 rho 32: -0.08477112361318034 rho 36: 0.11301403546356209
Aunque la serie es estacionaria, dado que las autocorrelaciones simples decrecen muy muy lento, hacemos una primera diferencia estacional.
DIFERENCIA ESTACIONAL¶
La diferencia estacional de orden $D$, con una periodicidad estacional $s$, se define de forma recursiva aplicando el operador $\nabla_s$ múltiples veces:
donde:
- $D$ es el orden de la diferencia estacional,
- $B$ es el operador de rezago tal que $B^k X_t = y_{X-k}$,
- $\nabla_s^D$ representa aplicar $D$ veces la diferencia estacional.
Para esta serie, el orden de la diferencia estacional es $D = 1$ y la periodicidad estacional es $s = 12$:
$$ \nabla_{12}^1 \text{T} (X_t) = \text{T} (X_t) - \text{T} (X_{t-12}) $$Sea entonces $Wt = \nabla_{12}^1 \text{T} (X_t) = \text{T} (X_t) - \text{T} (X_{t-12})$
PRIMERA DIFERENCIA ESTACIONAL¶
# Primera diferencia estacional con periodicidad 12
dzpre = zpre.diff(12).dropna()
# Visualizamos la serie
dzpre
195209
1954-02-01 -0.044262
1954-03-01 -0.621941
1954-04-01 0.357503
1954-05-01 1.805672
1954-06-01 3.170992
...
2008-08-01 0.762672
2008-09-01 0.396090
2008-10-01 0.048802
2008-11-01 0.039723
2008-12-01 0.000000
Length: 648, dtype: float64
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3)) # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) # Formato de fecha
plt.plot(dzpre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()
adfuller(dzpre)
(-10.871408688779647,
1.3655020494523785e-19,
11,
636,
{'1%': -3.4406737255613256,
'5%': -2.866095119842903,
'10%': -2.5691958123689727},
1448.5800337609196)
Estadístico ADF = -10.871
Valor-p = 1.37e-19
Número de rezagos = 11
Número de observaciones = 636
Valores críticos:
1% -> -3.4407
5% -> -2.8661
10% -> -2.5692
Log-likelihood = 9952.186
- El estadístico ADF es $-10.871$, menor que los valores críticos a los niveles del 1%.
- El valor-p es $1.37 × 10^{-19}$ mucho menor al nivel de significancia $\alpha$ = 0.05.
Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(dzpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(dzpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()
Se sigue notando la estacionalidad de periodicidad 12.
fac = FAC(len(dzpre), acf(dzpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r1: 0.0849542779953831 r5: -0.09455404757011154 r12: -0.414183329502961 r69: 0.13811523521698316 r73: -0.1125610007638793 r82: -0.09936326169869979
facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.08508558290727707 rho 5: -0.08389897753926695 rho 12: -0.435832204242704 rho 17: -0.12023508334631498 rho 24: -0.30097443805877616 rho 29: -0.10198920198664953 rho 36: -0.15950168604080855 rho 48: -0.14036369491687178 rho 57: -0.0916871776051421 rho 60: -0.12396066058923692 rho 61: 0.13800141862943816 rho 69: 0.10269095820413293 rho 72: -0.09462958643919762 rho 82: -0.08002594037040592 rho 85: -0.0819832551248085 rho 93: 0.0922037887338383 rho 109: 0.09582379018763905 rho 113: -0.09105599732819779 rho 119: 0.0917177196667104 rho 120: -0.1355872480803115
SEGUNDA DIFERENCIA ESTACIONAL¶
# 2 diferencia estacional con periodicidad 12
d2zpre = dzpre.diff(12).dropna()
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(d2zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(d2zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(d2zpre), acf(d2zpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r12: -0.6150182129084292 r57: -0.10922275236044279 r69: 0.16716305731426437 r73: -0.11365881167509073 r81: -0.12386243150076603 r82: -0.11601226324982332
facp = FACP(len(d2zpre), pacf(d2zpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 12: -0.6294838987246432 rho 17: -0.09052048792417981 rho 24: -0.5103011546739746 rho 29: -0.093217025773514 rho 36: -0.3551238626796082 rho 48: -0.2751156800940471 rho 57: -0.108916168798713 rho 60: -0.2379289221214959 rho 61: 0.20095345338550902 rho 65: -0.0834129475389058 rho 72: -0.29986531752257994 rho 73: 0.302496719075776 rho 74: -0.15214692323832477 rho 75: 0.12158013142041356 rho 77: -0.21501383907290988 rho 79: -0.08587576253968264 rho 80: 0.10575443412592936 rho 81: -0.11295475178675196 rho 83: 0.1412351900166441 rho 84: -0.47389037358187136 rho 85: 0.8052856123151159 rho 86: -2.7039482046544805 rho 87: -1.3736483987654273 rho 88: 1.0019992843716 rho 89: 235.89044017004235 rho 90: -0.9971937437125793 rho 91: 0.48021486078669995 rho 93: -0.20196306057247534 rho 94: 0.28532937400900626 rho 95: -0.23402172704695934 rho 101: 0.25176068174227584 rho 102: -0.5167154340590121 rho 103: 0.7690418805288879 rho 104: -1.3401560972346027 rho 105: -2.166656395549073 rho 106: 1.1988394283045958 rho 107: 2.0506797927028866 rho 108: -1.2599924429780232 rho 109: -2.603184507263249 rho 110: 0.5135974713732451 rho 111: 0.2590435314099863 rho 113: -0.3628723780659208 rho 114: 0.4207095680305757 rho 115: -0.3885886584070735 rho 116: 0.2697703099821219 rho 119: -0.12676720968125227 rho 120: 0.16418782499731185
TERCERA DIFERENCIA ESTACIONAL¶
# Primera diferencia estacional con periodicidad 12
d3zpre = d2zpre.diff(12).dropna()
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(d3zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(d3zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(d3zpre), acf(d3zpre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r12: -0.6989860750076662 r24: 0.19148304290921273 r57: -0.13077245780878802 r69: 0.18180783859402475 r81: -0.1446388376312966
facp = FACP(len(d3zpre), pacf(d3zpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 12: -0.7131085439918823 rho 24: -0.6165329897683706 rho 29: -0.08067623128823141 rho 36: -0.5171728973407944 rho 37: 0.1982055108755419 rho 40: 0.10015883076893528 rho 41: -0.09228473683010521 rho 47: 0.09394099253473005 rho 48: -0.4929991451837127 rho 49: 0.27281984306019774 rho 50: -0.09255185942814254 rho 51: -0.10940745619258534 rho 52: 0.34514660641128075 rho 53: -0.3969949656195508 rho 54: 0.33091157283957406 rho 55: -0.24725105044265083 rho 57: -0.1678397209757971 rho 58: -0.41066728720157253 rho 59: 1.4752062518993367 rho 60: 3.447013197287627 rho 61: -0.965945384029374 rho 62: 6.762979431106962 rho 63: 1.0827131685786548 rho 64: -0.8201163848967379 rho 65: 1.5200084948857941 rho 66: 1.3863162790341865 rho 67: -2.332119329595056 rho 68: -1.2962242829344859 rho 69: 0.5251991629964926 rho 70: 0.12535213223464997 rho 71: -0.43609447528833045 rho 72: 0.35806234464280917 rho 73: -0.1138999259876181 rho 74: -0.10235078255767406 rho 75: 0.31494875075818757 rho 76: -0.20483724449372892 rho 77: 0.0921074399984271 rho 78: -0.08928014769113522 rho 79: 0.0854845471954933 rho 80: -0.28927169831303046 rho 81: 0.20932339081585724 rho 82: -0.11854528891185322 rho 84: 0.19492618448545254 rho 87: 0.21562177268699964 rho 88: -0.19739536579003492 rho 89: 0.1375019109969488 rho 90: -0.21573215769368453 rho 91: 0.13649794688884304 rho 92: -0.4130856229291036 rho 93: 0.5769402061332968 rho 94: -1.1748169850119188 rho 95: -6.0784018693922 rho 96: 0.9057557611685417 rho 97: -0.5526621327578615 rho 98: 0.3652101944200118 rho 99: -0.14934082325742695 rho 100: -0.20792999327433961 rho 101: 0.32248566533747575 rho 102: -0.43340585213601085 rho 103: 0.35823263199853145 rho 104: -0.16127048470991623 rho 105: -0.09034092543015515 rho 106: 0.8633154494572077 rho 107: -9.435548190300198 rho 108: -1.0889219702203907 rho 109: 1.4801001973243877 rho 110: 2.1038789559878754 rho 111: -0.7905244792702399 rho 112: 0.6723557185188808 rho 113: -0.9378271953920853 rho 114: 5.770582312925257 rho 115: 0.947663800814546 rho 116: 3.6394621713184363 rho 117: -0.8146958649632807 rho 118: 0.18972318689108425 rho 119: 0.4988023348063774 rho 120: -0.48455031039030716
COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶
# Desviación estándar de la serie
stds = [np.std(zpre), np.std(dzpre), np.std(d2zpre), np.std(d3zpre)]
# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()
Nos quedamos con la primera diferencia estacional, porque tiene una mejor forma su FAC y FACP, además de tener menor varianza.
SERIE CON TRANSFORMACIÓN¶
Definición matemática de la transformación Yeo-Johnson, como se encuentra en la literatura estadística (propuesta por Ingram Olkin y I. Paul Yeo & R. J. Johnson en 2000).
Sea $x \in \mathbb{R}$ un valor (puede ser negativo, cero o positivo). La transformación $T(x; \lambda)$ se define como:
$$ T(x; \lambda) = \begin{cases} \frac{[(x + 1)^\lambda - 1]}{\lambda} & \text{si } x \geq 0, \, \lambda \ne 0 \\ \log(x + 1) & \text{si } x \geq 0, \, \lambda = 0 \\ -\frac{[(-x + 1)^{2 - \lambda} - 1]}{2 - \lambda} & \text{si } x < 0, \, \lambda \ne 2 \\ -\log(-x + 1) & \text{si } x < 0, \, \lambda = 2 \end{cases} $$ESTIMADOR PUNTUAL¶
from sklearn.preprocessing import PowerTransformer
import numpy as np
X = pre.values.reshape(-1, 1)
# Ajuste original para obtener lambda estimado
pt = PowerTransformer(method='yeo-johnson', standardize=True)
ypre = pt.fit_transform(X)
lambda_est = pt.lambdas_[0]
print("Lambda estimado:", pt.lambdas_) # Obtener el valor lambda estimado
Lambda estimado: [0.23849441]
INTERVALO DE CONFIANZA¶
# Bootstrap para el IC de lambda
n_boot = 1000
lambdas_boot = []
for _ in range(n_boot):
sample = np.random.choice(X.flatten(), size=len(X), replace=True).reshape(-1, 1)
pt_boot = PowerTransformer(method='yeo-johnson', standardize=False)
pt_boot.fit(sample)
lambdas_boot.append(pt_boot.lambdas_[0])
# Calcular percentiles para el intervalo de confianza al 95%
ci_lower = np.percentile(lambdas_boot, 2.5)
ci_upper = np.percentile(lambdas_boot, 97.5)
print(f"Estimación de λ: {lambda_est:.4f}")
print(f"IC 95% para λ: ({ci_lower:.4f}, {ci_upper:.4f})")
Estimación de λ: 0.2385 IC 95% para λ: (0.2097, 0.2693)
# Gráfico opcional del histograma de los λ bootstrap
plt.figure(figsize=(10, 6))
plt.fill_betweenx([0, 120], ci_lower, ci_upper, color=sns.color_palette()[4], alpha=0.15)
plt.hist(lambdas_boot, bins=30, alpha=0.7, edgecolor='white')
plt.axvline(ci_lower, label=f'2.5% = {ci_lower:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(ci_upper, label=f'97.5% = {ci_upper:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(lambda_est, label=f'λ original = {lambda_est:.4f}', color=sns.color_palette()[1], linewidth=3)
plt.title('Distribución Bootstrap de λ')
plt.xlabel('λ')
plt.ylabel('Frecuencia')
plt.legend()
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(0.01))
plt.savefig('imagenes/02-lambda-yj-bootstrapping.svg', bbox_inches='tight')
plt.show()
SERIE TRANSFORMADA¶
plt.figure()
plt.plot(date[:-12], ypre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Yeo-Johnson Estandarizado Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-precipitacion_yeo_johnson.svg', bbox_inches='tight')
plt.show()
ESTACIONARIEDAD¶
from statsmodels.tsa.stattools import adfuller
adfuller(ypre)
(-6.885194656270816,
1.4004074124320326e-09,
16,
643,
{'1%': -3.440560883168159,
'5%': -2.8660454146233434,
'10%': -2.569169329058723},
1139.5328800342036)
Estadístico ADF = -6.885
Valor-p = 1.40e-09
Número de rezagos = 16
Número de observaciones = 643
Valores críticos:
1% -> -3.4406
5% -> -2.8660
10% -> -2.5692
Log-likelihood = 3742.248
- El estadístico ADF es $-6.885$, menor que los valores críticos a los niveles del 1%.
- El valor-p es $1.40 × 10^{-9}$ mucho menor al nivel de significancia $\alpha$ = 0.05.
Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.
AUTOCORRELACIONES¶
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import * # Este es el archivo que se usa para el curso
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
axs[0].set_xlabel('Retraso')
axs[0].set_ylabel('ACF')
plot_pacf(ypre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].set_xlabel('Retraso')
axs[1].set_ylabel('PACF')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Transformada de la Precipitación Mensual en Puebla', fontsize=20, weight='bold')
plt.savefig('imagenes/04-acf_pacf.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(ypre), acf(ypre, nlags=48)[1:] )
Valores de autocorrelacion significativos: r1: 0.6331507478065659 r2: 0.3325589848089995 r4: -0.3429619382624403 r5: -0.6128788441981868 r6: -0.7155687723661789 r7: -0.5667143985229659 r8: -0.28172958031763784 r10: 0.34639100385471927 r11: 0.6139734232526982 r12: 0.6833003454277807 r13: 0.5238197525884694 r14: 0.23526370655336934 r16: -0.37665249496039976 r17: -0.6133232329123165 r18: -0.6522050881633413 r19: -0.49864844118107393 r22: 0.3581357207867556 r23: 0.6051041828734448 r24: 0.632952748472827 r25: 0.4787907060076014 r28: -0.3830626619493817 r29: -0.5771240646474324 r30: -0.6219569734501011 r31: -0.44116167635231773 r34: 0.3905997714697265 r35: 0.5842211455257492 r36: 0.606315712870404 r37: 0.41306326727063364 r40: -0.39053448733579516 r41: -0.5677702589671872 r42: -0.5689341833781766 r43: -0.3981104753267791 r46: 0.3759283566696746 r47: 0.5704749093905213 r48: 0.5490665087242427
facp = FACP(len(ypre), pacf(ypre, nlags=48)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.6341115228411734 rho 2: -0.11461335623878138 rho 3: -0.3000454427793004 rho 4: -0.33135817233932896 rho 5: -0.38109982701434264 rho 6: -0.32479269111279224 rho 7: -0.08061163300684875 rho 11: 0.2000309079008965 rho 12: 0.1532581642408276 rho 17: -0.12292322538615731 rho 22: -0.08358889406877636 rho 23: 0.12630684698989092 rho 30: -0.1028445923709159 rho 47: 0.09892532485023141
DIFERENCIAS ESTACIONALES¶
Por la FAC y la FACP, la serie transformada también necesita diferencias estacionales.
PRIMERA DIFERENCIA ESTACIONAL¶
# Primera diferencia estacional con periodicidad 12
ypre_series = pd.Series(ypre.flatten(), index=pre.index) # Convertimos a serie de tiempo
dypre = ypre_series.diff(12).dropna() # Diferencia estacional
# Visualizamos la serie
dypre
195209
1954-02-01 -0.109427
1954-03-01 -1.152261
1954-04-01 1.611746
1954-05-01 1.819524
1954-06-01 1.900692
...
2008-08-01 0.606731
2008-09-01 0.156192
2008-10-01 0.052883
2008-11-01 0.312361
2008-12-01 0.000000
Length: 648, dtype: float64
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3)) # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) # Formato de fecha
plt.plot(dypre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Transformación Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()
adfuller(dypre)
(-9.030402747697376,
5.464639337632109e-15,
16,
631,
{'1%': -3.440755866431696,
'5%': -2.86613130039063,
'10%': -2.569215089800357},
1308.464428619442)
Estadístico ADF = -10.871
Valor-p = 1.37e-19
Número de rezagos = 11
Número de observaciones = 636
Valores críticos:
1% -> -3.4407
5% -> -2.8661
10% -> -2.5692
Log-likelihood = 9952.186
- El estadístico ADF es $-10.871$, menor que los valores críticos a los niveles del 1%.
- El valor-p es $1.37 × 10^{-19}$ mucho menor al nivel de significancia $\alpha$ = 0.05.
Por lo tanto, se rechaza la hipótesis nula, es decir, se concluye que la serie es estacionaria.
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(dypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(dypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r1: 0.1881206416574384 r2: 0.12601374062799692 r6: -0.10023577930193152 r12: -0.40245032162240263 r33: 0.09802046981427195 r45: -0.11985655705785987
facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.18841139999075743 rho 2: 0.09425081653205467 rho 6: -0.0848567418000625 rho 12: -0.4319892618507997 rho 13: 0.08136260875048439 rho 17: -0.11385209887045368 rho 22: -0.08401445717662205 rho 24: -0.2809725408252081 rho 25: 0.10755279636113212 rho 30: -0.09673228306484968 rho 33: 0.1105360450418002 rho 36: -0.14602967562386657 rho 42: -0.11769525648448687 rho 46: -0.09378255882994088 rho 48: -0.16710968828081296 rho 50: 0.09151867173223789 rho 55: 0.09278816335947565 rho 60: -0.09971581325023068 rho 72: -0.13002571917786468 rho 78: -0.07964448430687596 rho 84: -0.07986957194068202 rho 92: -0.13966165907823416 rho 93: 0.11600081821072342 rho 94: -0.1221857229746919 rho 113: -0.08611614522897519 rho 118: -0.08199830520693245 rho 120: -0.11627799482364788
SEGUNDA DIFERENCIA ESTACIONAL¶
# 2 diferencia estacional con periodicidad 12
d2ypre = dypre.diff(12).dropna()
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(d2ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(d2ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(d2ypre), acf(d2ypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r12: -0.615018212908429 r57: -0.10922275236044277 r69: 0.16716305731426437 r73: -0.1136588116750907 r81: -0.12386243150076602 r82: -0.11601226324982332
facp = FACP(len(d2ypre), pacf(d2ypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 12: -0.6294838987246429 rho 17: -0.09052048792417974 rho 24: -0.510301154673974 rho 29: -0.09321702577351394 rho 36: -0.35512386267960755 rho 48: -0.2751156800940455 rho 57: -0.10891616879871288 rho 60: -0.23792892212149455 rho 61: 0.20095345338550874 rho 65: -0.08341294753890602 rho 72: -0.29986531752257806 rho 73: 0.30249671907577647 rho 74: -0.1521469232383268 rho 75: 0.12158013142041249 rho 77: -0.21501383907291427 rho 79: -0.08587576253968643 rho 80: 0.10575443412592932 rho 81: -0.11295475178675143 rho 83: 0.1412351900166407 rho 84: -0.4738903735818693 rho 85: 0.8052856123151381 rho 86: -2.7039482046551293 rho 87: -1.373648398765204 rho 88: 1.0019992843716292 rho 89: 235.89044016524818 rho 90: -0.9971937437125115 rho 91: 0.4802148607866808 rho 93: -0.20196306057246763 rho 94: 0.2853293740090054 rho 95: -0.23402172704695937 rho 101: 0.25176068174227473 rho 102: -0.5167154340590162 rho 103: 0.769041880528836 rho 104: -1.3401560972344024 rho 105: -2.166656395548571 rho 106: 1.198839428304655 rho 107: 2.0506797927018146 rho 108: -1.259992442977829 rho 109: -2.6031845072646265 rho 110: 0.5135974713735127 rho 111: 0.2590435314096445 rho 113: -0.36287237806607014 rho 114: 0.4207095680306412 rho 115: -0.3885886584071046 rho 116: 0.26977030998213086 rho 119: -0.12676720968126257 rho 120: 0.16418782499731296
TERCERA DIFERENCIA ESTACIONAL¶
# Primera diferencia estacional con periodicidad 12
d3ypre = d2ypre.diff(12).dropna()
fig, axs = plt.subplots(1, 2, figsize=(16, 4))
plot_acf(d3ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))
# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
plot_pacf(d3ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))
plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')
plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()
fac = FAC(len(d3ypre), acf(d3ypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r1: 0.14142525988694732 r11: -0.08618353748437811 r12: -0.7006965212651104 r13: -0.12299455598817151 r24: 0.21088909322401359 r33: 0.14799260300636627 r34: 0.124164682223292 r45: -0.13840321103959075
facp = FACP(len(d3ypre), pacf(d3ypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.14165226672464704 rho 11: -0.08572854737384034 rho 12: -0.7145741931403762 rho 14: 0.10244937563798918 rho 21: -0.0974412193826169 rho 22: -0.11958849854177123 rho 24: -0.590999784962588 rho 25: 0.1155406109341616 rho 36: -0.42596863520035666 rho 37: 0.18734873402939595 rho 42: -0.1733829621796275 rho 46: -0.1132553928924406 rho 48: -0.39708277000220765 rho 49: 0.17397913808655488 rho 54: -0.3269257232330364 rho 55: 0.21559931119807146 rho 56: -0.19500988553521115 rho 58: -0.22348004365902227 rho 59: 0.25219826387781036 rho 60: -0.7375993663060407 rho 61: 2.2537091363236095 rho 62: 1.650176227081085 rho 63: -0.6968349029252933 rho 64: 0.3967431321320999 rho 65: -0.2123951680489383 rho 66: -0.20497281303308482 rho 68: 0.23621005007258603 rho 69: -0.8189484098245967 rho 70: 1.3011193720883465 rho 71: 1.4316135793533165 rho 72: -5.898535785106428 rho 73: -1.3199340565173563 rho 74: -0.6772780595533806 rho 75: 0.5245892189329673 rho 76: -0.2815596779482279 rho 77: -0.4521464433566569 rho 78: 0.30058897281376856 rho 79: -0.4164535343317919 rho 80: -0.11859502105379839 rho 81: -0.15366212304252277 rho 82: -0.2922622187005957 rho 83: -0.4928118125453261 rho 84: -0.14977995949488318 rho 85: -0.8927516420039372 rho 86: -5.928900128261838 rho 87: 1.2117900757453872 rho 88: 0.5145451003205691 rho 89: 0.24638896698250676 rho 90: 0.5283271371105461 rho 92: 0.3314944963217833 rho 93: 0.11506584057595354 rho 94: 0.21220949431027086 rho 95: -0.12891904223431166 rho 96: 0.5141635439947064 rho 97: -0.6188892690196736 rho 98: 1.7101159305305729 rho 99: 2.12107036942252 rho 100: -0.6062698236020745 rho 101: 0.46103150897488143 rho 102: -0.1569944090737856 rho 103: 0.2535125462534724 rho 104: -0.21007258043287413 rho 105: 0.3130566565661735 rho 107: 0.1210008208606222 rho 109: 0.1042542157962568 rho 110: -0.17221115272098989 rho 111: 0.2646586117130243 rho 112: -0.11554768465356624 rho 114: 0.18496589221568 rho 115: 0.10819324030267996 rho 116: -0.15857098929802355 rho 117: 0.3781749655329778 rho 118: -0.1243216624045143 rho 120: 0.20596573821321837
COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶
# Desviación estándar de la serie
stds = [np.std(ypre), np.std(dypre), np.std(d2ypre), np.std(d3ypre)]
# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()
Nos quedamos con la primera diferencia estacional, porque tiene una mejor forma su FAC y FACP, además de tener menor varianza.
MODELADO¶
Lo que se muestra a continuación fue realizado para la serie sin transformar, solamente estandarizada.
En 02-----Encontrar-los-Mejores-Modelos se utilizá la libreria de pmdarima para intentar diferentes combinaciones de modelos y encontrar sus $\text{AIC}$.
Se obtuvieron los siguientes resultados:
(1) ARIMA(1,0,0)(0,1,1)[12], AIC=1430.397 (2) ARIMA(0,0,1)(0,1,1)[12], AIC=1431.266 (3) ARIMA(1,0,1)(0,1,1)[12], AIC=1432.133 (4) ARIMA(2,0,0)(0,1,1)[12], AIC=1432.157 (5) ARIMA(1,0,0)(0,1,2)[12], AIC=1432.368 (6) ARIMA(1,0,0)(1,1,1)[12], AIC=1432.373 (7) ARIMA(0,0,2)(0,1,1)[12], AIC=1432.664 (8) ARIMA(0,0,1)(0,1,2)[12], AIC=1433.240 (9) ARIMA(0,0,1)(1,1,1)[12], AIC=1433.245 (10) ARIMA(2,0,1)(0,1,1)[12], AIC=1434.127 (11) ARIMA(0,0,0)(0,1,1)[12], AIC=1445.334 (12) ARIMA(1,0,0)(1,1,0)[12], AIC=1523.141 (13) ARIMA(0,0,1)(1,1,0)[12], AIC=1523.219 (14) ARIMA(0,0,1)(0,1,0)[12], AIC=1652.849 (15) ARIMA(1,0,0)(0,1,0)[12], AIC=1652.984 (16) ARIMA(0,0,0)(0,1,0)[12], AIC=1655.671
Ningun modelo cumplió los supuestos. 02-01-Encontrar-Modelos-que-Cumplan
Por la FAC y FACP, se creyó que podía ser un $\text{ARIMA}(1, 0, 1) \times (0, 1, 1, 12)$, pero este modelo, con la prueba de Ljung-Box, tenía dependencia entre los residuos. Ver 03-----Encontrar-Modelos-que-Cumplan, 1 MODELO.
Se propuso aumentar el número de parámetros, se pusieron parámetros "intuitivos", es decir, la precipitación depende de los últimos 5 meses, de los últimos 4 años, y hay residuos en los últimos dos meses, llegando al primer modelo que cumplía la independecia en los RESIDUOS siendo el $\text{ARIMA}(5, 0, 2) \times (4, 1, 0, 12)$. Sin embargo, no era estacionario ni invertible. Ver 03-----Encontrar-Modelos-que-Cumplan, 2 MODELO.
Se fue variando el número de parámetros de los modelos, notando curiosamente, que para esta serie, que los modelos que eran estacionarios eran los que tenían su parte del polinomio autorregresivo completa y no tenían parte autorregresiva estacional (y viceversa), lo mismo para la parte media móvil y la invertibilidad. Tomando en cuenta esta observación, se hicieron las notebooks de 03-01, ..., 03-04, donde se varían, por ejemplo, los $p$, manteniendo $P=0$, lo mismo para la parte media móvil. Lo malo, es que para los modelos tuvieran independencia en los residuos, los polinomios tenían que ser bastante grandes, lo que podía hacer que no se cumpliera el principio de parsimonia.
En las notebooks 03-01, ..., 03-04 se proponen multiples modelos, que, dada la Observación (2), que dice que para esta serie los modelos que cumplen esa forma son admisibles, se siguió el siguiente procedimiento:
- Se fijan los valores de $\text{ARIMA}(p, 0, q) \times (0, 1, 0, 12)$ donde $p = 0, 1, 2, \dots$ y $q = 0, 1, 2, \dots$.
- Se fijan los valores de $\text{ARIMA}(0, 0, q) \times (P, 1, 0, 12)$ donde $P = 0, 1, 2, \dots$ y $q = 0, 1, 2, \dots$.
- Se fijan los valores de $\text{ARIMA}(p, 0, 0) \times (0, 1, Q, 12)$ donde $p = 0, 1, 2, \dots$ y $Q = 0, 1, 2, \dots$.
- Se fijan los valores de $\text{ARIMA}(0, 0, 0) \times (P, 1, Q, 12)$ donde $P = 0, 1, 2, \dots$ y $Q = 0, 1, 2, \dots$.
RESULTADOS de 03-01-Encontrar-Modelo-pq:
(1) ARIMA(0,0,15)(0,1,0)[12] AIC=1370.81 (2) ARIMA(0,0,12)(0,1,0)[12] AIC=1370.92 (3) ARIMA(3,0,12)(0,1,0)[12] AIC=1372.25 (4) ARIMA(0,0,16)(0,1,0)[12] AIC=1372.52 (5) ARIMA(1,0,15)(0,1,0)[12] AIC=1372.60 (6) ARIMA(1,0,12)(0,1,0)[12] AIC=1372.69 (7) ARIMA(0,0,13)(0,1,0)[12] AIC=1372.78 (8) ARIMA(3,0,13)(0,1,0)[12] AIC=1372.87 (9) ARIMA(4,0,12)(0,1,0)[12] AIC=1373.41 (10) ARIMA(1,0,16)(0,1,0)[12] AIC=1373.59 (11) ARIMA(1,0,13)(0,1,0)[12] AIC=1373.76 (12) ARIMA(3,0,14)(0,1,0)[12] AIC=1374.01 (13) ARIMA(4,0,13)(0,1,0)[12] AIC=1374.12 (14) ARIMA(5,0,12)(0,1,0)[12] AIC=1374.34 (15) ARIMA(0,0,14)(0,1,0)[12] AIC=1374.44 (16) ARIMA(2,0,12)(0,1,0)[12] AIC=1374.49 (17) ARIMA(2,0,12)(0,1,0)[12] AIC=1374.49 (18) ARIMA(0,0,17)(0,1,0)[12] AIC=1374.49 (19) ARIMA(2,0,15)(0,1,0)[12] AIC=1374.62 (20) ARIMA(2,0,15)(0,1,0)[12] AIC=1374.62 (21) ARIMA(1,0,17)(0,1,0)[12] AIC=1374.97 (22) ARIMA(4,0,14)(0,1,0)[12] AIC=1375.08 (23) ARIMA(2,0,13)(0,1,0)[12] AIC=1375.16 (24) ARIMA(2,0,13)(0,1,0)[12] AIC=1375.16 (25) ARIMA(1,0,14)(0,1,0)[12] AIC=1375.26 (26) ARIMA(5,0,13)(0,1,0)[12] AIC=1375.43 (27) ARIMA(6,0,12)(0,1,0)[12] AIC=1375.94 (28) ARIMA(6,0,12)(0,1,0)[12] AIC=1375.94 (29) ARIMA(2,0,16)(0,1,0)[12] AIC=1375.97 (30) ARIMA(2,0,16)(0,1,0)[12] AIC=1375.97 (31) ARIMA(3,0,15)(0,1,0)[12] AIC=1376.37 (32) ARIMA(7,0,12)(0,1,0)[12] AIC=1376.46 (33) ARIMA(2,0,14)(0,1,0)[12] AIC=1376.63 (34) ARIMA(2,0,14)(0,1,0)[12] AIC=1376.63 (35) ARIMA(5,0,14)(0,1,0)[12] AIC=1377.04 (36) ARIMA(6,0,13)(0,1,0)[12] AIC=1377.52 (37) ARIMA(6,0,13)(0,1,0)[12] AIC=1377.52 (38) ARIMA(4,0,15)(0,1,0)[12] AIC=1377.61 (39) ARIMA(2,0,17)(0,1,0)[12] AIC=1377.67 (40) ARIMA(2,0,17)(0,1,0)[12] AIC=1377.67 (41) ARIMA(3,0,16)(0,1,0)[12] AIC=1377.71 (42) ARIMA(7,0,13)(0,1,0)[12] AIC=1378.15 (43) ARIMA(8,0,12)(0,1,0)[12] AIC=1378.47 (44) ARIMA(8,0,12)(0,1,0)[12] AIC=1378.47 (45) ARIMA(9,0,12)(0,1,0)[12] AIC=1378.79 (46) ARIMA(6,0,14)(0,1,0)[12] AIC=1378.97 (47) ARIMA(6,0,14)(0,1,0)[12] AIC=1378.97 (48) ARIMA(4,0,16)(0,1,0)[12] AIC=1379.46 (49) ARIMA(3,0,17)(0,1,0)[12] AIC=1379.51 (50) ARIMA(5,0,15)(0,1,0)[12] AIC=1379.68 (51) ARIMA(8,0,13)(0,1,0)[12] AIC=1380.15 (52) ARIMA(8,0,13)(0,1,0)[12] AIC=1380.15 (53) ARIMA(7,0,14)(0,1,0)[12] AIC=1380.53 (54) ARIMA(10,0,12)(0,1,0)[12] AIC=1380.62 (55) ARIMA(5,0,16)(0,1,0)[12] AIC=1380.90 (56) ARIMA(4,0,17)(0,1,0)[12] AIC=1381.04 (57) ARIMA(6,0,15)(0,1,0)[12] AIC=1381.49 (58) ARIMA(6,0,15)(0,1,0)[12] AIC=1381.49 (59) ARIMA(10,0,13)(0,1,0)[12] AIC=1381.66 (60) ARIMA(9,0,13)(0,1,0)[12] AIC=1381.75 (61) ARIMA(8,0,14)(0,1,0)[12] AIC=1382.31 (62) ARIMA(8,0,14)(0,1,0)[12] AIC=1382.31 (63) ARIMA(5,0,17)(0,1,0)[12] AIC=1382.35 (64) ARIMA(11,0,12)(0,1,0)[12] AIC=1382.48 (65) ARIMA(9,0,14)(0,1,0)[12] AIC=1382.58 (66) ARIMA(6,0,16)(0,1,0)[12] AIC=1382.64 (67) ARIMA(6,0,16)(0,1,0)[12] AIC=1382.64 (68) ARIMA(7,0,15)(0,1,0)[12] AIC=1382.73 (69) ARIMA(10,0,14)(0,1,0)[12] AIC=1383.53 (70) ARIMA(11,0,13)(0,1,0)[12] AIC=1383.82 (71) ARIMA(6,0,17)(0,1,0)[12] AIC=1384.31 (72) ARIMA(6,0,17)(0,1,0)[12] AIC=1384.31 (73) ARIMA(7,0,16)(0,1,0)[12] AIC=1384.33 (74) ARIMA(8,0,15)(0,1,0)[12] AIC=1384.52 (75) ARIMA(8,0,15)(0,1,0)[12] AIC=1384.52 (76) ARIMA(9,0,15)(0,1,0)[12] AIC=1384.53 (77) ARIMA(13,0,12)(0,1,0)[12] AIC=1384.56 (78) ARIMA(11,0,14)(0,1,0)[12] AIC=1384.73 (79) ARIMA(13,0,13)(0,1,0)[12] AIC=1385.12 (80) ARIMA(10,0,15)(0,1,0)[12] AIC=1385.45 (81) ARIMA(12,0,13)(0,1,0)[12] AIC=1385.52 (82) ARIMA(8,0,16)(0,1,0)[12] AIC=1385.88 (83) ARIMA(8,0,16)(0,1,0)[12] AIC=1385.88 (84) ARIMA(12,0,12)(0,1,0)[12] AIC=1386.15 (85) ARIMA(7,0,17)(0,1,0)[12] AIC=1386.37 (86) ARIMA(11,0,15)(0,1,0)[12] AIC=1386.42 (87) ARIMA(9,0,16)(0,1,0)[12] AIC=1386.43 (88) ARIMA(13,0,14)(0,1,0)[12] AIC=1386.89 (89) ARIMA(12,0,14)(0,1,0)[12] AIC=1386.90 (90) ARIMA(10,0,16)(0,1,0)[12] AIC=1387.11 (91) ARIMA(14,0,13)(0,1,0)[12] AIC=1387.36 (92) ARIMA(14,0,12)(0,1,0)[12] AIC=1387.85 (93) ARIMA(11,0,16)(0,1,0)[12] AIC=1387.97 (94) ARIMA(8,0,17)(0,1,0)[12] AIC=1388.02 (95) ARIMA(8,0,17)(0,1,0)[12] AIC=1388.02 (96) ARIMA(15,0,12)(0,1,0)[12] AIC=1388.42 (97) ARIMA(13,0,15)(0,1,0)[12] AIC=1388.76 (98) ARIMA(12,0,15)(0,1,0)[12] AIC=1388.84 (99) ARIMA(14,0,14)(0,1,0)[12] AIC=1388.89 (100) ARIMA(9,0,17)(0,1,0)[12] AIC=1388.95 (101) ARIMA(10,0,17)(0,1,0)[12] AIC=1389.04 (102) ARIMA(15,0,13)(0,1,0)[12] AIC=1389.30 (103) ARIMA(17,0,12)(0,1,0)[12] AIC=1389.57 (104) ARIMA(11,0,17)(0,1,0)[12] AIC=1389.72 (105) ARIMA(16,0,12)(0,1,0)[12] AIC=1389.87 (106) ARIMA(12,0,16)(0,1,0)[12] AIC=1390.33 (107) ARIMA(15,0,14)(0,1,0)[12] AIC=1390.73 (108) ARIMA(14,0,15)(0,1,0)[12] AIC=1390.78 (109) ARIMA(17,0,13)(0,1,0)[12] AIC=1390.81 (110) ARIMA(16,0,13)(0,1,0)[12] AIC=1390.83 (111) ARIMA(13,0,16)(0,1,0)[12] AIC=1391.95 (112) ARIMA(12,0,17)(0,1,0)[12] AIC=1392.26 (113) ARIMA(13,0,17)(0,1,0)[12] AIC=1392.28 (114) ARIMA(17,0,14)(0,1,0)[12] AIC=1392.49 (115) ARIMA(16,0,14)(0,1,0)[12] AIC=1392.69 (116) ARIMA(15,0,15)(0,1,0)[12] AIC=1392.85 (117) ARIMA(14,0,16)(0,1,0)[12] AIC=1392.97 (118) ARIMA(16,0,15)(0,1,0)[12] AIC=1394.72 (119) ARIMA(17,0,15)(0,1,0)[12] AIC=1394.84 (120) ARIMA(15,0,16)(0,1,0)[12] AIC=1395.10 (121) ARIMA(14,0,17)(0,1,0)[12] AIC=1395.38 (122) ARIMA(15,0,17)(0,1,0)[12] AIC=1396.22 (123) ARIMA(17,0,16)(0,1,0)[12] AIC=1396.75 (124) ARIMA(16,0,16)(0,1,0)[12] AIC=1397.15 (125) ARIMA(16,0,17)(0,1,0)[12] AIC=1398.28 (126) ARIMA(17,0,17)(0,1,0)[12] AIC=1399.30 (127) ARIMA(9,0,11)(0,1,0)[12] AIC=1405.22 (128) ARIMA(11,0,10)(0,1,0)[12] AIC=1405.63 (129) ARIMA(11,0,11)(0,1,0)[12] AIC=1406.36 (130) ARIMA(10,0,11)(0,1,0)[12] AIC=1411.83 (131) ARIMA(14,0,11)(0,1,0)[12] AIC=1413.61 (132) ARIMA(13,0,11)(0,1,0)[12] AIC=1415.15 (133) ARIMA(8,0,11)(0,1,0)[12] AIC=1415.83 (134) ARIMA(8,0,11)(0,1,0)[12] AIC=1415.83 (135) ARIMA(15,0,11)(0,1,0)[12] AIC=1418.52 (136) ARIMA(16,0,11)(0,1,0)[12] AIC=1418.96 (137) ARIMA(12,0,11)(0,1,0)[12] AIC=1419.07 (138) ARIMA(7,0,11)(0,1,0)[12] AIC=1420.16 (139) ARIMA(15,0,10)(0,1,0)[12] AIC=1427.03 (140) ARIMA(16,0,10)(0,1,0)[12] AIC=1428.10 (141) ARIMA(16,0,10)(0,1,0)[12] AIC=1428.10 (142) ARIMA(14,0,10)(0,1,0)[12] AIC=1429.08 (143) ARIMA(4,0,11)(0,1,0)[12] AIC=1430.74 (144) ARIMA(3,0,11)(0,1,0)[12] AIC=1435.07 (145) ARIMA(17,0,10)(0,1,0)[12] AIC=1436.99 (146) ARIMA(12,0,10)(0,1,0)[12] AIC=1437.18 (147) ARIMA(5,0,11)(0,1,0)[12] AIC=1443.88 (148) ARIMA(13,0,10)(0,1,0)[12] AIC=1444.13 (149) ARIMA(6,0,11)(0,1,0)[12] AIC=1449.08 (150) ARIMA(6,0,11)(0,1,0)[12] AIC=1449.08 (151) ARIMA(16,0,9)(0,1,0)[12] AIC=1450.39 (152) ARIMA(17,0,9)(0,1,0)[12] AIC=1450.46 (153) ARIMA(2,0,11)(0,1,0)[12] AIC=1451.97 (154) ARIMA(2,0,11)(0,1,0)[12] AIC=1451.97 (155) ARIMA(16,0,8)(0,1,0)[12] AIC=1460.41 (156) ARIMA(17,0,8)(0,1,0)[12] AIC=1464.61 (157) ARIMA(16,0,7)(0,1,0)[12] AIC=1469.14 (158) ARIMA(17,0,7)(0,1,0)[12] AIC=1469.22 (159) ARIMA(17,0,6)(0,1,0)[12] AIC=1471.02 (160) ARIMA(9,0,10)(0,1,0)[12] AIC=1471.14 (161) ARIMA(8,0,10)(0,1,0)[12] AIC=1478.49 (162) ARIMA(8,0,10)(0,1,0)[12] AIC=1478.49 (163) ARIMA(16,0,5)(0,1,0)[12] AIC=1479.55 (164) ARIMA(16,0,6)(0,1,0)[12] AIC=1480.98 (165) ARIMA(10,0,10)(0,1,0)[12] AIC=1495.41 (166) ARIMA(17,0,5)(0,1,0)[12] AIC=1496.61 (167) ARIMA(6,0,10)(0,1,0)[12] AIC=1497.53 (168) ARIMA(6,0,10)(0,1,0)[12] AIC=1497.53 (169) ARIMA(1,0,11)(0,1,0)[12] AIC=1503.08 (170) ARIMA(5,0,10)(0,1,0)[12] AIC=1526.92 (171) ARIMA(10,0,6)(0,1,0)[12] AIC=1541.36 (172) ARIMA(4,0,10)(0,1,0)[12] AIC=1542.28 (173) ARIMA(7,0,10)(0,1,0)[12] AIC=1546.44 (174) ARIMA(3,0,10)(0,1,0)[12] AIC=1552.06 (175) ARIMA(2,0,10)(0,1,0)[12] AIC=1552.72 (176) ARIMA(2,0,10)(0,1,0)[12] AIC=1552.72 (177) ARIMA(2,0,10)(0,1,0)[12] AIC=1552.72 (178) ARIMA(2,0,9)(0,1,0)[12] AIC=1555.97 (179) ARIMA(2,0,8)(0,1,0)[12] AIC=1559.09 (180) ARIMA(1,0,8)(0,1,0)[12] AIC=1593.01 (181) ARIMA(0,0,11)(0,1,0)[12] AIC=1599.58 (182) ARIMA(2,0,5)(0,1,0)[12] AIC=1601.11 (183) ARIMA(2,0,6)(0,1,0)[12] AIC=1621.74 (184) ARIMA(0,0,9)(0,1,0)[12] AIC=1621.79 (185) ARIMA(0,0,10)(0,1,0)[12] AIC=1622.73 (186) ARIMA(0,0,10)(0,1,0)[12] AIC=1622.73 (187) ARIMA(1,0,9)(0,1,0)[12] AIC=1623.45 (188) ARIMA(1,0,10)(0,1,0)[12] AIC=1624.20 (189) ARIMA(1,0,10)(0,1,0)[12] AIC=1624.20 (190) ARIMA(1,0,5)(0,1,0)[12] AIC=1627.56 (191) ARIMA(1,0,6)(0,1,0)[12] AIC=1630.49 (192) ARIMA(0,0,7)(0,1,0)[12] AIC=1630.59 (193) ARIMA(0,0,8)(0,1,0)[12] AIC=1630.74 (194) ARIMA(1,0,7)(0,1,0)[12] AIC=1631.92 (195) ARIMA(2,0,7)(0,1,0)[12] AIC=1633.08 (196) ARIMA(0,0,6)(0,1,0)[12] AIC=1635.97 (197) ARIMA(0,0,5)(0,1,0)[12] AIC=1651.91
RESULTADOS de 03-02-Encontrar-Modelo-Pq:
(1) ARIMA(0,0,15)(0,1,0)[12] AIC=1370.81 (2) ARIMA(0,0,12)(0,1,0)[12] AIC=1370.92 (3) ARIMA(0,0,12)(1,1,0)[12] AIC=1372.79 (4) ARIMA(0,0,15)(1,1,0)[12] AIC=1373.17 (5) ARIMA(0,0,12)(2,1,0)[12] AIC=1373.37 (6) ARIMA(0,0,13)(1,1,0)[12] AIC=1373.86 (7) ARIMA(0,0,12)(4,1,0)[12] AIC=1375.93 (8) ARIMA(0,0,14)(2,1,0)[12] AIC=1377.49 (9) ARIMA(0,0,11)(3,1,0)[12] AIC=1427.35 (10) ARIMA(0,0,10)(3,1,0)[12] AIC=1430.51 (11) ARIMA(0,0,11)(2,1,0)[12] AIC=1443.12 (12) ARIMA(0,0,10)(2,1,0)[12] AIC=1446.97 (13) ARIMA(0,0,11)(1,1,0)[12] AIC=1511.21 (14) ARIMA(0,0,6)(1,1,0)[12] AIC=1514.47 (15) ARIMA(0,0,8)(1,1,0)[12] AIC=1516.15 (16) ARIMA(0,0,10)(1,1,0)[12] AIC=1516.98
RESULTADOS de 03-03-Encontrar-Modelo-pQ:
(1) ARIMA(12,0,0)(0,1,4)[12] AIC=1367.04 (2) ARIMA(13,0,0)(0,1,4)[12] AIC=1367.37 (3) ARIMA(12,0,0)(0,1,3)[12] AIC=1368.24 (4) ARIMA(12,0,0)(0,1,5)[12] AIC=1369.25 (5) ARIMA(15,0,0)(0,1,4)[12] AIC=1371.49 (6) ARIMA(12,0,0)(0,1,1)[12] AIC=1382.66 (7) ARIMA(15,0,0)(0,1,1)[12] AIC=1387.68 (8) ARIMA(6,0,0)(0,1,1)[12] AIC=1391.61 (9) ARIMA(11,0,0)(0,1,4)[12] AIC=1392.63 (10) ARIMA(10,0,0)(0,1,4)[12] AIC=1398.90 (11) ARIMA(3,0,0)(0,1,1)[12] AIC=1431.92 (12) ARIMA(0,0,0)(0,1,1)[12] AIC=1443.59 (13) ARIMA(12,0,0)(0,1,0)[12] AIC=1522.91 (14) ARIMA(6,0,0)(0,1,0)[12] AIC=1650.81
RESULTADOS de 03-04-Encontrar-Modelo-PQ:
(1) ARIMA(0,0,0)(0,1,3)[12] AIC=1443.50 (2) ARIMA(0,0,0)(0,1,1)[12] AIC=1443.59 (3) ARIMA(0,0,0)(2,1,1)[12] AIC=1444.11 (4) ARIMA(0,0,0)(1,1,2)[12] AIC=1445.15 (5) ARIMA(0,0,0)(1,1,3)[12] AIC=1445.28 (6) ARIMA(0,0,0)(0,1,4)[12] AIC=1445.49 (7) ARIMA(0,0,0)(0,1,2)[12] AIC=1445.59 (8) ARIMA(0,0,0)(1,1,1)[12] AIC=1445.59 (9) ARIMA(0,0,0)(4,1,3)[12] AIC=1445.70 (10) ARIMA(0,0,0)(4,1,4)[12] AIC=1445.74 (11) ARIMA(0,0,0)(2,1,2)[12] AIC=1445.92 (12) ARIMA(0,0,0)(4,1,2)[12] AIC=1445.95 (13) ARIMA(0,0,0)(3,1,1)[12] AIC=1446.07 (14) ARIMA(0,0,0)(2,1,3)[12] AIC=1446.38 (15) ARIMA(0,0,0)(1,1,4)[12] AIC=1446.57 (16) ARIMA(0,0,0)(4,1,1)[12] AIC=1446.70 (17) ARIMA(0,0,0)(3,1,2)[12] AIC=1446.99 (18) ARIMA(0,0,0)(3,1,3)[12] AIC=1447.78 (19) ARIMA(0,0,0)(2,1,4)[12] AIC=1447.80 (20) ARIMA(0,0,0)(4,1,0)[12] AIC=1448.24 (21) ARIMA(0,0,0)(3,1,4)[12] AIC=1449.21 (22) ARIMA(0,0,0)(3,1,0)[12] AIC=1454.13 (23) ARIMA(0,0,0)(2,1,0)[12] AIC=1464.14 (24) ARIMA(0,0,0)(1,1,0)[12] AIC=1526.64 (25) ARIMA(0,0,0)(0,1,0)[12] AIC=1653.67
Estos son todos los modelos propuestos, ordenados por su AIC:
(1) ARIMA(12,0,0)(0,1,4)[12], AIC=1367.04 (2) ARIMA(13,0,0)(0,1,4)[12], AIC=1367.37 (3) ARIMA(12,0,0)(0,1,3)[12], AIC=1368.24 (4) ARIMA(12,0,0)(0,1,5)[12], AIC=1369.25 (5) ARIMA(0,0,15)(0,1,0)[12], AIC=1370.81 (6) ARIMA(0,0,15)(0,1,0)[12], AIC=1370.81 (7) ARIMA(0,0,12)(0,1,0)[12], AIC=1370.92 (8) ARIMA(0,0,12)(0,1,0)[12], AIC=1370.92 (9) ARIMA(15,0,0)(0,1,4)[12], AIC=1371.49 (10) ARIMA(3,0,12)(0,1,0)[12], AIC=1372.25 (11) ARIMA(0,0,16)(0,1,0)[12], AIC=1372.52 (12) ARIMA(1,0,15)(0,1,0)[12], AIC=1372.60 (13) ARIMA(1,0,12)(0,1,0)[12], AIC=1372.69 (14) ARIMA(0,0,13)(0,1,0)[12], AIC=1372.78 (15) ARIMA(0,0,12)(1,1,0)[12], AIC=1372.79 (16) ARIMA(3,0,13)(0,1,0)[12], AIC=1372.87 (17) ARIMA(0,0,15)(1,1,0)[12], AIC=1373.17 (18) ARIMA(0,0,12)(2,1,0)[12], AIC=1373.37 (19) ARIMA(4,0,12)(0,1,0)[12], AIC=1373.41 (20) ARIMA(1,0,16)(0,1,0)[12], AIC=1373.59 (21) ARIMA(1,0,13)(0,1,0)[12], AIC=1373.76 (22) ARIMA(0,0,13)(1,1,0)[12], AIC=1373.86 (23) ARIMA(3,0,14)(0,1,0)[12], AIC=1374.01 (24) ARIMA(4,0,13)(0,1,0)[12], AIC=1374.12 (25) ARIMA(5,0,12)(0,1,0)[12], AIC=1374.34 (26) ARIMA(0,0,14)(0,1,0)[12], AIC=1374.44 (27) ARIMA(2,0,12)(0,1,0)[12], AIC=1374.49 (28) ARIMA(2,0,12)(0,1,0)[12], AIC=1374.49 (29) ARIMA(0,0,17)(0,1,0)[12], AIC=1374.49 (30) ARIMA(2,0,15)(0,1,0)[12], AIC=1374.62 (31) ARIMA(2,0,15)(0,1,0)[12], AIC=1374.62 (32) ARIMA(1,0,17)(0,1,0)[12], AIC=1374.97 (33) ARIMA(4,0,14)(0,1,0)[12], AIC=1375.08 (34) ARIMA(2,0,13)(0,1,0)[12], AIC=1375.16 (35) ARIMA(2,0,13)(0,1,0)[12], AIC=1375.16 (36) ARIMA(1,0,14)(0,1,0)[12], AIC=1375.26 (37) ARIMA(5,0,13)(0,1,0)[12], AIC=1375.43 (38) ARIMA(0,0,12)(4,1,0)[12], AIC=1375.93 (39) ARIMA(6,0,12)(0,1,0)[12], AIC=1375.94 (40) ARIMA(6,0,12)(0,1,0)[12], AIC=1375.94 (41) ARIMA(2,0,16)(0,1,0)[12], AIC=1375.97 (42) ARIMA(2,0,16)(0,1,0)[12], AIC=1375.97 (43) ARIMA(3,0,15)(0,1,0)[12], AIC=1376.37 (44) ARIMA(7,0,12)(0,1,0)[12], AIC=1376.46 (45) ARIMA(2,0,14)(0,1,0)[12], AIC=1376.63 (46) ARIMA(2,0,14)(0,1,0)[12], AIC=1376.63 (47) ARIMA(5,0,14)(0,1,0)[12], AIC=1377.04 (48) ARIMA(0,0,14)(2,1,0)[12], AIC=1377.49 (49) ARIMA(6,0,13)(0,1,0)[12], AIC=1377.52 (50) ARIMA(6,0,13)(0,1,0)[12], AIC=1377.52 (51) ARIMA(4,0,15)(0,1,0)[12], AIC=1377.61 (52) ARIMA(2,0,17)(0,1,0)[12], AIC=1377.67 (53) ARIMA(2,0,17)(0,1,0)[12], AIC=1377.67 (54) ARIMA(3,0,16)(0,1,0)[12], AIC=1377.71 (55) ARIMA(7,0,13)(0,1,0)[12], AIC=1378.15 (56) ARIMA(8,0,12)(0,1,0)[12], AIC=1378.47 (57) ARIMA(8,0,12)(0,1,0)[12], AIC=1378.47 (58) ARIMA(9,0,12)(0,1,0)[12], AIC=1378.79 (59) ARIMA(6,0,14)(0,1,0)[12], AIC=1378.97 (60) ARIMA(6,0,14)(0,1,0)[12], AIC=1378.97 (61) ARIMA(4,0,16)(0,1,0)[12], AIC=1379.46 (62) ARIMA(3,0,17)(0,1,0)[12], AIC=1379.51 (63) ARIMA(5,0,15)(0,1,0)[12], AIC=1379.68 (64) ARIMA(8,0,13)(0,1,0)[12], AIC=1380.15 (65) ARIMA(8,0,13)(0,1,0)[12], AIC=1380.15 (66) ARIMA(7,0,14)(0,1,0)[12], AIC=1380.53 (67) ARIMA(10,0,12)(0,1,0)[12], AIC=1380.62 (68) ARIMA(5,0,16)(0,1,0)[12], AIC=1380.90 (69) ARIMA(4,0,17)(0,1,0)[12], AIC=1381.04 (70) ARIMA(6,0,15)(0,1,0)[12], AIC=1381.49 (71) ARIMA(6,0,15)(0,1,0)[12], AIC=1381.49 (72) ARIMA(10,0,13)(0,1,0)[12], AIC=1381.66 (73) ARIMA(9,0,13)(0,1,0)[12], AIC=1381.75 (74) ARIMA(8,0,14)(0,1,0)[12], AIC=1382.31 (75) ARIMA(8,0,14)(0,1,0)[12], AIC=1382.31 (76) ARIMA(5,0,17)(0,1,0)[12], AIC=1382.35 (77) ARIMA(11,0,12)(0,1,0)[12], AIC=1382.48 (78) ARIMA(9,0,14)(0,1,0)[12], AIC=1382.58 (79) ARIMA(6,0,16)(0,1,0)[12], AIC=1382.64 (80) ARIMA(6,0,16)(0,1,0)[12], AIC=1382.64 (81) ARIMA(12,0,0)(0,1,1)[12], AIC=1382.66 (82) ARIMA(7,0,15)(0,1,0)[12], AIC=1382.73 (83) ARIMA(10,0,14)(0,1,0)[12], AIC=1383.53 (84) ARIMA(11,0,13)(0,1,0)[12], AIC=1383.82 (85) ARIMA(6,0,17)(0,1,0)[12], AIC=1384.31 (86) ARIMA(6,0,17)(0,1,0)[12], AIC=1384.31 (87) ARIMA(7,0,16)(0,1,0)[12], AIC=1384.33 (88) ARIMA(8,0,15)(0,1,0)[12], AIC=1384.52 (89) ARIMA(8,0,15)(0,1,0)[12], AIC=1384.52 (90) ARIMA(9,0,15)(0,1,0)[12], AIC=1384.53 (91) ARIMA(13,0,12)(0,1,0)[12], AIC=1384.56 (92) ARIMA(11,0,14)(0,1,0)[12], AIC=1384.73 (93) ARIMA(13,0,13)(0,1,0)[12], AIC=1385.12 (94) ARIMA(10,0,15)(0,1,0)[12], AIC=1385.45 (95) ARIMA(12,0,13)(0,1,0)[12], AIC=1385.52 (96) ARIMA(8,0,16)(0,1,0)[12], AIC=1385.88 (97) ARIMA(8,0,16)(0,1,0)[12], AIC=1385.88 (98) ARIMA(12,0,12)(0,1,0)[12], AIC=1386.15 (99) ARIMA(7,0,17)(0,1,0)[12], AIC=1386.37 (100) ARIMA(11,0,15)(0,1,0)[12], AIC=1386.42 (101) ARIMA(9,0,16)(0,1,0)[12], AIC=1386.43 (102) ARIMA(13,0,14)(0,1,0)[12], AIC=1386.89 (103) ARIMA(12,0,14)(0,1,0)[12], AIC=1386.90 (104) ARIMA(10,0,16)(0,1,0)[12], AIC=1387.11 (105) ARIMA(14,0,13)(0,1,0)[12], AIC=1387.36 (106) ARIMA(15,0,0)(0,1,1)[12], AIC=1387.68 (107) ARIMA(14,0,12)(0,1,0)[12], AIC=1387.85 (108) ARIMA(11,0,16)(0,1,0)[12], AIC=1387.97 (109) ARIMA(8,0,17)(0,1,0)[12], AIC=1388.02 (110) ARIMA(8,0,17)(0,1,0)[12], AIC=1388.02 (111) ARIMA(15,0,12)(0,1,0)[12], AIC=1388.42 (112) ARIMA(13,0,15)(0,1,0)[12], AIC=1388.76 (113) ARIMA(12,0,15)(0,1,0)[12], AIC=1388.84 (114) ARIMA(14,0,14)(0,1,0)[12], AIC=1388.89 (115) ARIMA(9,0,17)(0,1,0)[12], AIC=1388.95 (116) ARIMA(10,0,17)(0,1,0)[12], AIC=1389.04 (117) ARIMA(15,0,13)(0,1,0)[12], AIC=1389.30 (118) ARIMA(17,0,12)(0,1,0)[12], AIC=1389.57 (119) ARIMA(11,0,17)(0,1,0)[12], AIC=1389.72 (120) ARIMA(16,0,12)(0,1,0)[12], AIC=1389.87 (121) ARIMA(12,0,16)(0,1,0)[12], AIC=1390.33 (122) ARIMA(15,0,14)(0,1,0)[12], AIC=1390.73 (123) ARIMA(14,0,15)(0,1,0)[12], AIC=1390.78 (124) ARIMA(17,0,13)(0,1,0)[12], AIC=1390.81 (125) ARIMA(16,0,13)(0,1,0)[12], AIC=1390.83 (126) ARIMA(6,0,0)(0,1,1)[12], AIC=1391.61 (127) ARIMA(13,0,16)(0,1,0)[12], AIC=1391.95 (128) ARIMA(12,0,17)(0,1,0)[12], AIC=1392.26 (129) ARIMA(13,0,17)(0,1,0)[12], AIC=1392.28 (130) ARIMA(17,0,14)(0,1,0)[12], AIC=1392.49 (131) ARIMA(11,0,0)(0,1,4)[12], AIC=1392.63 (132) ARIMA(16,0,14)(0,1,0)[12], AIC=1392.69 (133) ARIMA(15,0,15)(0,1,0)[12], AIC=1392.85 (134) ARIMA(14,0,16)(0,1,0)[12], AIC=1392.97 (135) ARIMA(16,0,15)(0,1,0)[12], AIC=1394.72 (136) ARIMA(17,0,15)(0,1,0)[12], AIC=1394.84 (137) ARIMA(15,0,16)(0,1,0)[12], AIC=1395.10 (138) ARIMA(14,0,17)(0,1,0)[12], AIC=1395.38 (139) ARIMA(15,0,17)(0,1,0)[12], AIC=1396.22 (140) ARIMA(17,0,16)(0,1,0)[12], AIC=1396.75 (141) ARIMA(16,0,16)(0,1,0)[12], AIC=1397.15 (142) ARIMA(16,0,17)(0,1,0)[12], AIC=1398.28 (143) ARIMA(10,0,0)(0,1,4)[12], AIC=1398.90 (144) ARIMA(17,0,17)(0,1,0)[12], AIC=1399.30 (145) ARIMA(9,0,11)(0,1,0)[12], AIC=1405.22 (146) ARIMA(11,0,10)(0,1,0)[12], AIC=1405.63 (147) ARIMA(11,0,11)(0,1,0)[12], AIC=1406.36 (148) ARIMA(10,0,11)(0,1,0)[12], AIC=1411.83 (149) ARIMA(14,0,11)(0,1,0)[12], AIC=1413.61 (150) ARIMA(13,0,11)(0,1,0)[12], AIC=1415.15 (151) ARIMA(8,0,11)(0,1,0)[12], AIC=1415.83 (152) ARIMA(8,0,11)(0,1,0)[12], AIC=1415.83 (153) ARIMA(15,0,11)(0,1,0)[12], AIC=1418.52 (154) ARIMA(16,0,11)(0,1,0)[12], AIC=1418.96 (155) ARIMA(12,0,11)(0,1,0)[12], AIC=1419.07 (156) ARIMA(7,0,11)(0,1,0)[12], AIC=1420.16 (157) ARIMA(15,0,10)(0,1,0)[12], AIC=1427.03 (158) ARIMA(0,0,11)(3,1,0)[12], AIC=1427.35 (159) ARIMA(16,0,10)(0,1,0)[12], AIC=1428.10 (160) ARIMA(16,0,10)(0,1,0)[12], AIC=1428.10 (161) ARIMA(14,0,10)(0,1,0)[12], AIC=1429.08 (162) ARIMA(0,0,10)(3,1,0)[12], AIC=1430.51 (163) ARIMA(4,0,11)(0,1,0)[12], AIC=1430.74 (164) ARIMA(3,0,0)(0,1,1)[12], AIC=1431.92 (165) ARIMA(3,0,11)(0,1,0)[12], AIC=1435.07 (166) ARIMA(17,0,10)(0,1,0)[12], AIC=1436.99 (167) ARIMA(12,0,10)(0,1,0)[12], AIC=1437.18 (168) ARIMA(0,0,11)(2,1,0)[12], AIC=1443.12 (169) ARIMA(0,0,0)(0,1,3)[12], AIC=1443.50 (170) ARIMA(0,0,0)(0,1,1)[12], AIC=1443.59 (171) ARIMA(0,0,0)(0,1,1)[12], AIC=1443.59 (172) ARIMA(5,0,11)(0,1,0)[12], AIC=1443.88 (173) ARIMA(0,0,0)(2,1,1)[12], AIC=1444.11 (174) ARIMA(13,0,10)(0,1,0)[12], AIC=1444.13 (175) ARIMA(0,0,0)(1,1,2)[12], AIC=1445.15 (176) ARIMA(0,0,0)(1,1,3)[12], AIC=1445.28 (177) ARIMA(0,0,0)(0,1,4)[12], AIC=1445.49 (178) ARIMA(0,0,0)(0,1,2)[12], AIC=1445.59 (179) ARIMA(0,0,0)(1,1,1)[12], AIC=1445.59 (180) ARIMA(0,0,0)(4,1,3)[12], AIC=1445.70 (181) ARIMA(0,0,0)(4,1,4)[12], AIC=1445.74 (182) ARIMA(0,0,0)(2,1,2)[12], AIC=1445.92 (183) ARIMA(0,0,0)(4,1,2)[12], AIC=1445.95 (184) ARIMA(0,0,0)(3,1,1)[12], AIC=1446.07 (185) ARIMA(0,0,0)(2,1,3)[12], AIC=1446.38 (186) ARIMA(0,0,0)(1,1,4)[12], AIC=1446.57 (187) ARIMA(0,0,0)(4,1,1)[12], AIC=1446.70 (188) ARIMA(0,0,10)(2,1,0)[12], AIC=1446.97 (189) ARIMA(0,0,0)(3,1,2)[12], AIC=1446.99 (190) ARIMA(0,0,0)(3,1,3)[12], AIC=1447.78 (191) ARIMA(0,0,0)(2,1,4)[12], AIC=1447.80 (192) ARIMA(0,0,0)(4,1,0)[12], AIC=1448.24 (193) ARIMA(6,0,11)(0,1,0)[12], AIC=1449.08 (194) ARIMA(6,0,11)(0,1,0)[12], AIC=1449.08 (195) ARIMA(0,0,0)(3,1,4)[12], AIC=1449.21 (196) ARIMA(16,0,9)(0,1,0)[12], AIC=1450.39 (197) ARIMA(17,0,9)(0,1,0)[12], AIC=1450.46 (198) ARIMA(2,0,11)(0,1,0)[12], AIC=1451.97 (199) ARIMA(2,0,11)(0,1,0)[12], AIC=1451.97 (200) ARIMA(0,0,0)(3,1,0)[12], AIC=1454.13 (201) ARIMA(16,0,8)(0,1,0)[12], AIC=1460.41 (202) ARIMA(0,0,0)(2,1,0)[12], AIC=1464.14 (203) ARIMA(17,0,8)(0,1,0)[12], AIC=1464.61 (204) ARIMA(16,0,7)(0,1,0)[12], AIC=1469.14 (205) ARIMA(17,0,7)(0,1,0)[12], AIC=1469.22 (206) ARIMA(17,0,6)(0,1,0)[12], AIC=1471.02 (207) ARIMA(9,0,10)(0,1,0)[12], AIC=1471.14 (208) ARIMA(8,0,10)(0,1,0)[12], AIC=1478.49 (209) ARIMA(8,0,10)(0,1,0)[12], AIC=1478.49 (210) ARIMA(16,0,5)(0,1,0)[12], AIC=1479.55 (211) ARIMA(16,0,6)(0,1,0)[12], AIC=1480.98 (212) ARIMA(10,0,10)(0,1,0)[12], AIC=1495.41 (213) ARIMA(17,0,5)(0,1,0)[12], AIC=1496.61 (214) ARIMA(6,0,10)(0,1,0)[12], AIC=1497.53 (215) ARIMA(6,0,10)(0,1,0)[12], AIC=1497.53 (216) ARIMA(1,0,11)(0,1,0)[12], AIC=1503.08 (217) ARIMA(0,0,11)(1,1,0)[12], AIC=1511.21 (218) ARIMA(0,0,6)(1,1,0)[12], AIC=1514.47 (219) ARIMA(0,0,8)(1,1,0)[12], AIC=1516.15 (220) ARIMA(0,0,10)(1,1,0)[12], AIC=1516.98 (221) ARIMA(12,0,0)(0,1,0)[12], AIC=1522.91 (222) ARIMA(0,0,0)(1,1,0)[12], AIC=1526.64 (223) ARIMA(5,0,10)(0,1,0)[12], AIC=1526.92 (224) ARIMA(10,0,6)(0,1,0)[12], AIC=1541.36 (225) ARIMA(4,0,10)(0,1,0)[12], AIC=1542.28 (226) ARIMA(7,0,10)(0,1,0)[12], AIC=1546.44 (227) ARIMA(3,0,10)(0,1,0)[12], AIC=1552.06 (228) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72 (229) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72 (230) ARIMA(2,0,10)(0,1,0)[12], AIC=1552.72 (231) ARIMA(2,0,9)(0,1,0)[12], AIC=1555.97 (232) ARIMA(2,0,8)(0,1,0)[12], AIC=1559.09 (233) ARIMA(1,0,8)(0,1,0)[12], AIC=1593.01 (234) ARIMA(0,0,11)(0,1,0)[12], AIC=1599.58 (235) ARIMA(2,0,5)(0,1,0)[12], AIC=1601.11 (236) ARIMA(2,0,6)(0,1,0)[12], AIC=1621.74 (237) ARIMA(0,0,9)(0,1,0)[12], AIC=1621.79 (238) ARIMA(0,0,10)(0,1,0)[12], AIC=1622.73 (239) ARIMA(0,0,10)(0,1,0)[12], AIC=1622.73 (240) ARIMA(1,0,9)(0,1,0)[12], AIC=1623.45 (241) ARIMA(1,0,10)(0,1,0)[12], AIC=1624.20 (242) ARIMA(1,0,10)(0,1,0)[12], AIC=1624.20 (243) ARIMA(1,0,5)(0,1,0)[12], AIC=1627.56 (244) ARIMA(1,0,6)(0,1,0)[12], AIC=1630.49 (245) ARIMA(0,0,7)(0,1,0)[12], AIC=1630.59 (246) ARIMA(0,0,8)(0,1,0)[12], AIC=1630.74 (247) ARIMA(1,0,7)(0,1,0)[12], AIC=1631.92 (248) ARIMA(2,0,7)(0,1,0)[12], AIC=1633.08 (249) ARIMA(0,0,6)(0,1,0)[12], AIC=1635.97 (250) ARIMA(6,0,0)(0,1,0)[12], AIC=1650.81 (251) ARIMA(0,0,5)(0,1,0)[12], AIC=1651.91 (252) ARIMA(0,0,0)(0,1,0)[12], AIC=1653.67
MODELOS PROPUESTOS¶
from statsmodels.tsa.statespace.sarimax import SARIMAX
from verificacion_de_supuestos import * # Una funcion que hice para verificar si los parametros son estacionarios/invertibles
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import ttest_1samp
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.tools.tools import add_constant
from scipy.stats import jarque_bera
from statsmodels.stats.diagnostic import lilliefors
MODELOS PARA LA SERIE SIN TRANSFORMAR¶
1 MODELO¶
modelo=SARIMAX(zpre,
order=(12,0,0),
seasonal_order=(0,1,4,12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(12, 0, 0)x(0, 1, [1, 2, 3, 4], 12) | Log Likelihood | -666.522 |
| Date: | Sun, 27 Apr 2025 | AIC | 1367.044 |
| Time: | 04:48:12 | BIC | 1443.100 |
| Sample: | 0 | HQIC | 1396.549 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.1079 | 0.032 | 3.338 | 0.001 | 0.045 | 0.171 |
| ar.L2 | -0.0271 | 0.037 | -0.738 | 0.460 | -0.099 | 0.045 |
| ar.L3 | 0.0504 | 0.037 | 1.367 | 0.172 | -0.022 | 0.123 |
| ar.L4 | -0.0689 | 0.042 | -1.634 | 0.102 | -0.152 | 0.014 |
| ar.L5 | -0.1086 | 0.050 | -2.193 | 0.028 | -0.206 | -0.012 |
| ar.L6 | -0.1103 | 0.059 | -1.872 | 0.061 | -0.226 | 0.005 |
| ar.L7 | -0.0049 | 0.050 | -0.097 | 0.923 | -0.104 | 0.094 |
| ar.L8 | -0.0450 | 0.043 | -1.039 | 0.299 | -0.130 | 0.040 |
| ar.L9 | 0.0411 | 0.035 | 1.183 | 0.237 | -0.027 | 0.109 |
| ar.L10 | -0.0133 | 0.033 | -0.406 | 0.685 | -0.078 | 0.051 |
| ar.L11 | 0.0968 | 0.036 | 2.662 | 0.008 | 0.026 | 0.168 |
| ar.L12 | 0.5442 | 0.074 | 7.328 | 0.000 | 0.399 | 0.690 |
| ma.S.L12 | -1.3768 | 0.084 | -16.419 | 0.000 | -1.541 | -1.212 |
| ma.S.L24 | 0.3656 | 0.079 | 4.639 | 0.000 | 0.211 | 0.520 |
| ma.S.L36 | 0.1054 | 0.060 | 1.756 | 0.079 | -0.012 | 0.223 |
| ma.S.L48 | -0.0790 | 0.039 | -2.042 | 0.041 | -0.155 | -0.003 |
| sigma2 | 0.4378 | 0.024 | 18.022 | 0.000 | 0.390 | 0.485 |
| Ljung-Box (L1) (Q): | 0.34 | Jarque-Bera (JB): | 376.73 |
|---|---|---|---|
| Prob(Q): | 0.56 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.82 | Skew: | 0.89 |
| Prob(H) (two-sided): | 0.14 | Kurtosis: | 6.28 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
PRINCIPIO DE PARSIMONIA¶
parsimonia(modelo)
Hay coeficientes no significativos, no se cumple el principio de parsimonia
MODELO ADMISIBLE¶
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.10787169558218329, 'p2': 0.027085643770389415, 'p3': -0.050379964730704004, 'p4': 0.06894551834156623, 'p5': 0.10861700542171245, 'p6': 0.11029529000485752, 'p7': 0.00488096912398539, 'p8': 0.04497037748451452, 'p9': -0.04112819961848579, 'p10': 0.013332216884030874, 'p11': -0.09678291912170793, 'p12': -0.5442160802231802}
El grado del polinomio es: 12
Raíces del polinomio característico: [-1.10915274+0.j -0.92789955+0.52544644j -0.92789955-0.52544644j
-0.52132244+0.91355779j -0.52132244-0.91355779j -0.01548825+1.04638174j
-0.01548825-1.04638174j 0.53546751+0.91668064j 0.53546751-0.91668064j
1.05996383+0.j 0.86491762+0.50841197j 0.86491762-0.50841197j]
Módulo de las raíces: [1.10915274 1.06634494 1.06634494 1.05183883 1.05183883 1.04649636
1.04649636 1.06161615 1.06161615 1.05996383 1.00327724 1.00327724]
¿Las raíces están fuera del círculo unitario? True
El modelo es estacionario
:)
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.376803188196061, 'p24': 0.36563094250781014, 'p36': 0.10536384747112354, 'p48': -0.07904343456059729}
El grado del polinomio es: 48
Raíces del polinomio característico: [-1.05513205e+00+0.28272178j -1.05513205e+00-0.28272178j
-1.06072484e+00+0.06092667j -1.06072484e+00-0.06092667j
-1.00194881e+00+0.j -7.72410271e-01+0.77241027j
-7.72410271e-01-0.77241027j -9.49077989e-01+0.47759838j
-9.49077989e-01-0.47759838j -8.88151320e-01+0.58312646j
-8.88151320e-01-0.58312646j -8.67713125e-01+0.50097441j
-8.67713125e-01-0.50097441j -5.83126461e-01+0.88815132j
-5.83126461e-01-0.88815132j -4.77598375e-01+0.94907799j
-4.77598375e-01-0.94907799j -5.00974406e-01+0.86771313j
-5.00974406e-01-0.86771313j -2.82721781e-01+1.05513205j
-2.82721781e-01-1.05513205j -6.09266684e-02+1.06072484j
-6.09266684e-02-1.06072484j 6.09266684e-02+1.06072484j
6.09266684e-02-1.06072484j 4.53730117e-17+1.00194881j
4.53730117e-17-1.00194881j 2.82721781e-01+1.05513205j
2.82721781e-01-1.05513205j 4.77598375e-01+0.94907799j
4.77598375e-01-0.94907799j 5.83126461e-01+0.88815132j
5.83126461e-01-0.88815132j 5.00974406e-01+0.86771313j
5.00974406e-01-0.86771313j 7.72410271e-01+0.77241027j
7.72410271e-01-0.77241027j 1.05513205e+00+0.28272178j
1.05513205e+00-0.28272178j 1.06072484e+00+0.06092667j
1.06072484e+00-0.06092667j 1.00194881e+00+0.j
8.88151320e-01+0.58312646j 8.88151320e-01-0.58312646j
9.49077989e-01+0.47759838j 9.49077989e-01-0.47759838j
8.67713125e-01+0.50097441j 8.67713125e-01-0.50097441j]
Módulo de las raíces: [1.09235308 1.09235308 1.06247317 1.06247317 1.00194881 1.09235308
1.09235308 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
1.00194881 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
1.00194881 1.09235308 1.09235308 1.06247317 1.06247317 1.06247317
1.06247317 1.00194881 1.00194881 1.09235308 1.09235308 1.06247317
1.06247317 1.06247317 1.06247317 1.00194881 1.00194881 1.09235308
1.09235308 1.09235308 1.09235308 1.06247317 1.06247317 1.00194881
1.06247317 1.06247317 1.06247317 1.06247317 1.00194881 1.00194881]
¿Las raíces están fuera del círculo unitario? True
El modelo es invertible
:)
RESIDUOS INDEPENDIENTES¶
La prueba de Ljung-Box evalúa si los residuos de un modelo están autocorrelacionados, es decir, si existe dependencia temporal que el modelo no ha capturado. El resumen muestra lo siguiente:
Estadístico Q de Ljung-Box (L=1): 0.00
Valor-p asociado (Prob(Q)): 0.95
- El estadístico Q = 0.00 indica que no se detecta autocorrelación en los residuos al primer rezago.
- El valor-p = 0.95 es mucho mayor que el umbral común de significancia (α = 0.05).
Por lo tanto, no se rechaza la hipótesis nula de independencia, los residuos se comportan de manera independiente.
La prueba Ljung-Box evalúa la autocorrelación de los residuos, cuando estos residuos son residuales de un modelo ajustado, el estadístico debe corregirse para reflejar la pérdida de grados de libertad debido a la estimación de parámetros.
La prueba original de Ljung-Box asume que los datos son observaciones independientes sin parámetros estimados Cuando se usan residuos ajustados por un modelo con $k$ parámetros estimados, los grados de libertad del estadístico Q se reducen en función de estos parámetros. El argumento model_df en acorr_ljungbox() resta esos grados de libertad automáticamente.
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 1.410877 | NaN |
| 2 | 1.666736 | NaN |
| 3 | 1.668529 | NaN |
| 4 | 2.153928 | NaN |
| 5 | 3.756515 | NaN |
| 6 | 4.271537 | NaN |
| 7 | 4.271537 | NaN |
| 8 | 4.398558 | NaN |
| 9 | 4.673186 | NaN |
| 10 | 5.854609 | NaN |
| 11 | 6.235510 | NaN |
| 12 | 6.364705 | NaN |
| 13 | 8.908873 | NaN |
| 14 | 9.013624 | NaN |
| 15 | 9.115022 | NaN |
| 16 | 9.301973 | NaN |
| 17 | 9.622827 | 0.001922 |
| 18 | 9.634894 | 0.008087 |
| 19 | 9.648719 | 0.021801 |
| 20 | 9.650120 | 0.046752 |
| 21 | 10.177061 | 0.070371 |
| 22 | 10.403762 | 0.108646 |
| 23 | 10.403778 | 0.166823 |
| 24 | 10.406790 | 0.237627 |
| 25 | 10.406958 | 0.318555 |
| 26 | 10.956519 | 0.360916 |
| 27 | 11.577235 | 0.396244 |
| 28 | 11.768293 | 0.464463 |
| 29 | 12.219944 | 0.509696 |
| 30 | 12.223270 | 0.588378 |
| 31 | 12.229827 | 0.661556 |
| 32 | 14.130904 | 0.588962 |
| 33 | 14.189628 | 0.653636 |
| 34 | 16.160520 | 0.581346 |
| 35 | 16.186598 | 0.644794 |
| 36 | 16.663770 | 0.674686 |
| 37 | 16.699639 | 0.729145 |
| 38 | 17.046283 | 0.760802 |
| 39 | 17.392114 | 0.789577 |
| 40 | 18.027704 | 0.801662 |
| 41 | 19.922150 | 0.750897 |
| 42 | 20.004090 | 0.791363 |
| 43 | 20.038756 | 0.829126 |
| 44 | 20.096537 | 0.860920 |
| 45 | 20.910376 | 0.862394 |
| 46 | 21.002354 | 0.887814 |
| 47 | 21.052280 | 0.910525 |
| 48 | 21.957194 | 0.908534 |
Hay dependencia en los primeros dos residuos.
RESIDUOS CON MEDIA CERO¶
modelo.resid.mean()
-0.04082317526120878
ttest_1samp(modelo.resid, 0)
TtestResult(statistic=-1.5347957438291504, pvalue=0.12531394692788528, df=659)
La media de los residuos puede ser 0.
RESIDUOS CON VARIANZA CONSTANTE¶
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
(2.3754715860699527, 0.12325437239232831, 2.376827864684162, 0.12362813746633024)
Los residuos tienen varianza constante.
RESIDUOS CON DISTRIBUCIÓN NORMAL¶
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=417.01556109393744, pvalue=2.793968301284977e-91) (0.11546272377541411, 0.0009999999999998899)
Los residuos no siguen una distribución normal.
normal_relajacion(modelo.resid)
75.30% de los residuos están dentro de ±1σ (esperado ≈ 68%) 93.94% de los residuos están dentro de ±2σ (esperado ≈ 95%) 98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)
Pero son bastante parecidos.
GRÁFICO DE RESIDUOS¶
from matplotlib.collections import PathCollection
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)
axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
line.set_linewidth(0.5)
# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
line.set_alpha(0.2)
# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)
plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-01-analisis-de-residuos-m1.svg', bbox_inches='tight')
plt.show()
2 MODELO¶
modelo=SARIMAX(zpre,
order=(13,0,0),
seasonal_order=(0,1,4,12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) | Log Likelihood | -665.685 |
| Date: | Sun, 27 Apr 2025 | AIC | 1367.369 |
| Time: | 04:49:22 | BIC | 1447.899 |
| Sample: | 0 | HQIC | 1398.609 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.1369 | 0.035 | 3.963 | 0.000 | 0.069 | 0.205 |
| ar.L2 | -0.0190 | 0.036 | -0.522 | 0.602 | -0.090 | 0.052 |
| ar.L3 | 0.0501 | 0.037 | 1.363 | 0.173 | -0.022 | 0.122 |
| ar.L4 | -0.0697 | 0.042 | -1.662 | 0.097 | -0.152 | 0.013 |
| ar.L5 | -0.1117 | 0.050 | -2.254 | 0.024 | -0.209 | -0.015 |
| ar.L6 | -0.1101 | 0.060 | -1.839 | 0.066 | -0.227 | 0.007 |
| ar.L7 | -0.0104 | 0.051 | -0.205 | 0.838 | -0.110 | 0.089 |
| ar.L8 | -0.0511 | 0.044 | -1.170 | 0.242 | -0.137 | 0.034 |
| ar.L9 | 0.0331 | 0.035 | 0.944 | 0.345 | -0.036 | 0.102 |
| ar.L10 | -0.0093 | 0.033 | -0.281 | 0.778 | -0.074 | 0.055 |
| ar.L11 | 0.1034 | 0.037 | 2.802 | 0.005 | 0.031 | 0.176 |
| ar.L12 | 0.5385 | 0.080 | 6.699 | 0.000 | 0.381 | 0.696 |
| ar.L13 | -0.0460 | 0.041 | -1.115 | 0.265 | -0.127 | 0.035 |
| ma.S.L12 | -1.3756 | 0.100 | -13.725 | 0.000 | -1.572 | -1.179 |
| ma.S.L24 | 0.3674 | 0.086 | 4.275 | 0.000 | 0.199 | 0.536 |
| ma.S.L36 | 0.0973 | 0.060 | 1.635 | 0.102 | -0.019 | 0.214 |
| ma.S.L48 | -0.0784 | 0.039 | -1.994 | 0.046 | -0.155 | -0.001 |
| sigma2 | 0.4340 | 0.030 | 14.503 | 0.000 | 0.375 | 0.493 |
| Ljung-Box (L1) (Q): | 0.04 | Jarque-Bera (JB): | 347.08 |
|---|---|---|---|
| Prob(Q): | 0.85 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | 0.87 |
| Prob(H) (two-sided): | 0.17 | Kurtosis: | 6.14 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
PRINCIPIO DE PARSIMONIA¶
parsimonia(modelo)
Hay coeficientes no significativos, no se cumple el principio de parsimonia
MODELO ADMISIBLE¶
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.13691650540076578, 'p2': 0.01901377784903687, 'p3': -0.0500667165465603, 'p4': 0.069698229094628, 'p5': 0.11171395157927178, 'p6': 0.11009957522513744, 'p7': 0.010425775539985482, 'p8': 0.051137304909435626, 'p9': -0.03308758191153757, 'p10': 0.009267663653695928, 'p11': -0.1033592404899208, 'p12': -0.5385061478654415, 'p13': 0.046023546700947306}
El grado del polinomio es: 13
Raíces del polinomio característico: [11.88851581+0.j -1.09964507+0.j -0.92538665+0.52179453j
-0.92538665-0.52179453j -0.52533631+0.91031656j -0.52533631-0.91031656j
-0.02218345+1.04735764j -0.02218345-1.04735764j 0.5303036 +0.92193348j
0.5303036 -0.92193348j 1.06751362+0.j 0.86474363+0.50810694j
0.86474363-0.50810694j]
Módulo de las raíces: [11.88851581 1.09964507 1.06236057 1.06236057 1.05102544 1.05102544
1.04759254 1.04759254 1.06357099 1.06357099 1.06751362 1.00297268
1.00297268]
¿Las raíces están fuera del círculo unitario? True
El modelo es estacionario
:)
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.3756005193776482, 'p24': 0.3674412910433979, 'p36': 0.09734997881875754, 'p48': -0.07835471066657075}
El grado del polinomio es: 48
Raíces del polinomio característico: [-1.05628158e+00+0.2830298j -1.05628158e+00-0.2830298j
-1.06074970e+00+0.06258621j -1.06074970e+00-0.06258621j
-1.00135998e+00+0.j -7.73251787e-01+0.77325179j
-7.73251787e-01-0.77325179j -9.49929287e-01+0.4761736j
-9.49929287e-01-0.4761736j -8.87343081e-01+0.58457609j
-8.87343081e-01-0.58457609j -8.67203179e-01+0.50067999j
-8.67203179e-01-0.50067999j -5.84576092e-01+0.88734308j
-5.84576092e-01-0.88734308j -4.76173604e-01+0.94992929j
-4.76173604e-01-0.94992929j -5.00679989e-01+0.86720318j
-5.00679989e-01-0.86720318j -2.83029797e-01+1.05628158j
-2.83029797e-01-1.05628158j -6.25862060e-02+1.0607497j
-6.25862060e-02-1.0607497j 6.25862060e-02+1.0607497j
6.25862060e-02-1.0607497j 2.16421979e-16+1.00135998j
2.16421979e-16-1.00135998j 2.83029797e-01+1.05628158j
2.83029797e-01-1.05628158j 4.76173604e-01+0.94992929j
4.76173604e-01-0.94992929j 5.84576092e-01+0.88734308j
5.84576092e-01-0.88734308j 5.00679989e-01+0.86720318j
5.00679989e-01-0.86720318j 7.73251787e-01+0.77325179j
7.73251787e-01-0.77325179j 1.05628158e+00+0.2830298j
1.05628158e+00-0.2830298j 1.06074970e+00+0.06258621j
1.06074970e+00-0.06258621j 1.00135998e+00+0.j
8.87343081e-01+0.58457609j 8.87343081e-01-0.58457609j
9.49929287e-01+0.4761736j 9.49929287e-01-0.4761736j
8.67203179e-01+0.50067999j 8.67203179e-01-0.50067999j]
Módulo de las raíces: [1.09354316 1.09354316 1.06259444 1.06259444 1.00135998 1.09354316
1.09354316 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
1.00135998 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
1.00135998 1.09354316 1.09354316 1.06259444 1.06259444 1.06259444
1.06259444 1.00135998 1.00135998 1.09354316 1.09354316 1.06259444
1.06259444 1.06259444 1.06259444 1.00135998 1.00135998 1.09354316
1.09354316 1.09354316 1.09354316 1.06259444 1.06259444 1.00135998
1.06259444 1.06259444 1.06259444 1.06259444 1.00135998 1.00135998]
¿Las raíces están fuera del círculo unitario? True
El modelo es invertible
:)
RESIDUOS INDEPENDIENTES¶
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 0.169555 | NaN |
| 2 | 0.204491 | NaN |
| 3 | 0.207502 | NaN |
| 4 | 0.727124 | NaN |
| 5 | 2.068440 | NaN |
| 6 | 2.452464 | NaN |
| 7 | 2.514097 | NaN |
| 8 | 2.809591 | NaN |
| 9 | 3.353825 | NaN |
| 10 | 4.344480 | NaN |
| 11 | 4.491955 | NaN |
| 12 | 4.756059 | NaN |
| 13 | 5.816217 | NaN |
| 14 | 6.013503 | NaN |
| 15 | 6.142261 | NaN |
| 16 | 6.454918 | NaN |
| 17 | 6.774386 | NaN |
| 18 | 6.783618 | 0.009200 |
| 19 | 6.786944 | 0.033592 |
| 20 | 6.809256 | 0.078232 |
| 21 | 7.196463 | 0.125863 |
| 22 | 7.428032 | 0.190703 |
| 23 | 7.441854 | 0.281909 |
| 24 | 7.447962 | 0.383773 |
| 25 | 7.506681 | 0.483077 |
| 26 | 8.263074 | 0.507868 |
| 27 | 8.898238 | 0.541788 |
| 28 | 9.010612 | 0.620913 |
| 29 | 9.533962 | 0.656772 |
| 30 | 9.544856 | 0.730662 |
| 31 | 9.545139 | 0.794626 |
| 32 | 11.177851 | 0.739888 |
| 33 | 11.194537 | 0.797321 |
| 34 | 13.243411 | 0.719745 |
| 35 | 13.251323 | 0.776431 |
| 36 | 13.698894 | 0.800953 |
| 37 | 13.804935 | 0.840247 |
| 38 | 14.263347 | 0.858024 |
| 39 | 14.567144 | 0.880135 |
| 40 | 15.189601 | 0.887710 |
| 41 | 16.919636 | 0.852069 |
| 42 | 16.949758 | 0.883608 |
| 43 | 16.977019 | 0.909776 |
| 44 | 17.083844 | 0.929034 |
| 45 | 17.970345 | 0.926894 |
| 46 | 18.039168 | 0.943470 |
| 47 | 18.078955 | 0.957241 |
| 48 | 18.964369 | 0.955485 |
Hay dependencia en los primeros dos residuos.
RESIDUOS CON MEDIA CERO¶
modelo.resid.mean()
-0.04384260119451572
ttest_1samp(modelo.resid, 0)
TtestResult(statistic=-1.6518512196908124, pvalue=0.09904125403907901, df=659)
La media de los residuos puede ser 0.
RESIDUOS CON VARIANZA CONSTANTE¶
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
(2.2724202057175913, 0.1316941016400431, 2.273361405689974, 0.13209411797812018)
Los residuos tienen varianza constante.
RESIDUOS CON DISTRIBUCIÓN NORMAL¶
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=380.3670597932282, pvalue=2.537128226453911e-83) (0.11315000343530657, 0.0009999999999998899)
Los residuos no siguen una distribución normal.
normal_relajacion(modelo.resid)
75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%) 93.79% de los residuos están dentro de ±2σ (esperado ≈ 95%) 98.94% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)
Pero son bastante parecidos.
GRÁFICO DE RESIDUOS¶
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)
axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
line.set_linewidth(0.5)
# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas
# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
# line.set_markerfacecolor('green')
# line.set_markeredgecolor('white')
line.set_alpha(0.2)
# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)
plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(13, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-02-analisis-de-residuos-m2.svg', bbox_inches='tight')
plt.show()
3 MODELO¶
modelo=SARIMAX(zpre,
order=(12,0,0),
seasonal_order=(0,1,2,12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(12, 0, 0)x(0, 1, [1, 2], 12) | Log Likelihood | -668.499 |
| Date: | Sun, 27 Apr 2025 | AIC | 1366.998 |
| Time: | 04:49:32 | BIC | 1434.107 |
| Sample: | 0 | HQIC | 1393.032 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.1241 | 0.033 | 3.797 | 0.000 | 0.060 | 0.188 |
| ar.L2 | -0.0278 | 0.037 | -0.746 | 0.455 | -0.101 | 0.045 |
| ar.L3 | 0.0593 | 0.036 | 1.629 | 0.103 | -0.012 | 0.131 |
| ar.L4 | -0.0777 | 0.043 | -1.819 | 0.069 | -0.161 | 0.006 |
| ar.L5 | -0.1077 | 0.050 | -2.137 | 0.033 | -0.206 | -0.009 |
| ar.L6 | -0.1154 | 0.060 | -1.914 | 0.056 | -0.234 | 0.003 |
| ar.L7 | -0.0107 | 0.050 | -0.214 | 0.831 | -0.108 | 0.087 |
| ar.L8 | -0.0498 | 0.044 | -1.129 | 0.259 | -0.136 | 0.037 |
| ar.L9 | 0.0450 | 0.036 | 1.263 | 0.207 | -0.025 | 0.115 |
| ar.L10 | -0.0087 | 0.033 | -0.263 | 0.793 | -0.074 | 0.056 |
| ar.L11 | 0.1075 | 0.035 | 3.042 | 0.002 | 0.038 | 0.177 |
| ar.L12 | 0.5090 | 0.073 | 6.931 | 0.000 | 0.365 | 0.653 |
| ma.S.L12 | -1.3374 | 0.082 | -16.355 | 0.000 | -1.498 | -1.177 |
| ma.S.L24 | 0.3550 | 0.075 | 4.760 | 0.000 | 0.209 | 0.501 |
| sigma2 | 0.4394 | 0.022 | 20.157 | 0.000 | 0.397 | 0.482 |
| Ljung-Box (L1) (Q): | 0.11 | Jarque-Bera (JB): | 385.54 |
|---|---|---|---|
| Prob(Q): | 0.74 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.82 | Skew: | 0.91 |
| Prob(H) (two-sided): | 0.15 | Kurtosis: | 6.31 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
PRINCIPIO DE PARSIMONIA¶
parsimonia(modelo)
Hay coeficientes no significativos, no se cumple el principio de parsimonia
MODELO ADMISIBLE¶
estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.12409609720857327, 'p2': 0.02784753160299208, 'p3': -0.05931370060464292, 'p4': 0.07773595850485018, 'p5': 0.10769845417883445, 'p6': 0.1153846887487895, 'p7': 0.010651434003279378, 'p8': 0.04983749921571501, 'p9': -0.04496683215549213, 'p10': 0.008718941011892536, 'p11': -0.10747007955137028, 'p12': -0.5090305660410539}
El grado del polinomio es: 12
Raíces del polinomio característico: [-1.12406749+0.j -0.93496549+0.52952748j -0.93496549-0.52952748j
-0.52508333+0.91815394j -0.52508333-0.91815394j -0.01640444+1.05341344j
-0.01640444-1.05341344j 0.53677332+0.92317355j 0.53677332-0.92317355j
1.06348079+0.j 0.8644098 +0.50791775j 0.8644098 -0.50791775j]
Módulo de las raíces: [1.12406749 1.07450445 1.07450445 1.05769521 1.05769521 1.05354116
1.05354116 1.06788342 1.06788342 1.06348079 1.00258902 1.00258902]
¿Las raíces están fuera del círculo unitario? True
El modelo es estacionario
:)
invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -1.3373896067612656, 'p24': 0.3549819568174022}
El grado del polinomio es: 24
Raíces del polinomio característico: [-1.08759159e+00+0.j -1.00234440e+00+0.j
-9.41881945e-01+0.54379579j -9.41881945e-01-0.54379579j
-8.68055716e-01+0.5011722j -8.68055716e-01-0.5011722j
-5.43795795e-01+0.94188195j -5.43795795e-01-0.94188195j
-5.01172202e-01+0.86805572j -5.01172202e-01-0.86805572j
-7.34340403e-16+1.08759159j -7.34340403e-16-1.08759159j
9.15927823e-17+1.0023444j 9.15927823e-17-1.0023444j
5.43795795e-01+0.94188195j 5.43795795e-01-0.94188195j
5.01172202e-01+0.86805572j 5.01172202e-01-0.86805572j
9.41881945e-01+0.54379579j 9.41881945e-01-0.54379579j
1.08759159e+00+0.j 1.00234440e+00+0.j
8.68055716e-01+0.5011722j 8.68055716e-01-0.5011722j ]
Módulo de las raíces: [1.08759159 1.0023444 1.08759159 1.08759159 1.0023444 1.0023444
1.08759159 1.08759159 1.0023444 1.0023444 1.08759159 1.08759159
1.0023444 1.0023444 1.08759159 1.08759159 1.0023444 1.0023444
1.08759159 1.08759159 1.08759159 1.0023444 1.0023444 1.0023444 ]
¿Las raíces están fuera del círculo unitario? True
El modelo es invertible
:)
RESIDUOS INDEPENDIENTES¶
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 0.852590 | NaN |
| 2 | 1.026959 | NaN |
| 3 | 1.042580 | NaN |
| 4 | 1.330274 | NaN |
| 5 | 2.702490 | NaN |
| 6 | 2.993605 | NaN |
| 7 | 3.102314 | NaN |
| 8 | 3.441180 | NaN |
| 9 | 3.616131 | NaN |
| 10 | 4.475840 | NaN |
| 11 | 4.662731 | NaN |
| 12 | 4.990946 | NaN |
| 13 | 7.916099 | NaN |
| 14 | 8.223545 | NaN |
| 15 | 8.469831 | 0.003611 |
| 16 | 8.751457 | 0.012579 |
| 17 | 9.104454 | 0.027934 |
| 18 | 9.135493 | 0.057801 |
| 19 | 9.135514 | 0.103780 |
| 20 | 9.139624 | 0.165877 |
| 21 | 9.745116 | 0.203483 |
| 22 | 10.168911 | 0.253370 |
| 23 | 10.181542 | 0.335989 |
| 24 | 10.580242 | 0.391141 |
| 25 | 10.587480 | 0.478441 |
| 26 | 11.470408 | 0.489092 |
| 27 | 12.051651 | 0.523413 |
| 28 | 12.201265 | 0.590143 |
| 29 | 12.738996 | 0.622450 |
| 30 | 12.744188 | 0.691362 |
| 31 | 12.749628 | 0.752779 |
| 32 | 14.953637 | 0.665148 |
| 33 | 14.972683 | 0.724332 |
| 34 | 17.346330 | 0.630384 |
| 35 | 17.349312 | 0.689726 |
| 36 | 18.361216 | 0.684370 |
| 37 | 18.361346 | 0.737605 |
| 38 | 18.844803 | 0.760215 |
| 39 | 19.320376 | 0.781506 |
| 40 | 19.893036 | 0.796598 |
| 41 | 21.815050 | 0.746715 |
| 42 | 21.912901 | 0.785308 |
| 43 | 21.985648 | 0.820783 |
| 44 | 22.126466 | 0.849404 |
| 45 | 23.141840 | 0.843941 |
| 46 | 23.250141 | 0.870261 |
| 47 | 23.337109 | 0.893652 |
| 48 | 23.478694 | 0.912247 |
Hay dependencia en los primeros dos residuos.
RESIDUOS CON MEDIA CERO¶
modelo.resid.mean()
-0.03975299918077584
ttest_1samp(modelo.resid, 0)
TtestResult(statistic=-1.4907248613412258, pvalue=0.1365122225435534, df=659)
La media de los residuos puede ser 0.
RESIDUOS CON VARIANZA CONSTANTE¶
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
(2.3477675031968714, 0.12546257027135715, 2.3490090062319307, 0.1258433113034253)
Los residuos tienen varianza constante.
RESIDUOS CON DISTRIBUCIÓN NORMAL¶
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=421.16509310160836, pvalue=3.508827740565011e-92) (0.10731720679113188, 0.0009999999999998899)
Los residuos no siguen una distribución normal.
normal_relajacion(modelo.resid)
75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%) 94.09% de los residuos están dentro de ±2σ (esperado ≈ 95%) 98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)
Pero son bastante parecidos.
GRÁFICO DE RESIDUOS¶
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)
axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
line.set_linewidth(0.5)
# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas
# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
# line.set_markerfacecolor('green')
# line.set_markeredgecolor('white')
line.set_alpha(0.2)
# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)
plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
OTROS MODELOS QUE SE CREYÓ QUE PODÍAN SER¶
Luego de ver la FAC y FACP:
fac = FAC(len(dzpre), acf(dzpre, nlags=12*12)[1:] )
Valores de autocorrelacion significativos: r1: 0.08495427799538309 r5: -0.09455404757011152 r12: -0.41418332950296094 r69: 0.13811523521698316 r73: -0.11256100076387932 r82: -0.09936326169869977 r139: 0.10081286249593209 r144: -0.103722152043645
facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.08508558290727704 rho 5: -0.08389897753926691 rho 12: -0.43583220424270375 rho 17: -0.12023508334631484 rho 24: -0.30097443805877616 rho 29: -0.10198920198664936 rho 36: -0.15950168604080825 rho 48: -0.14036369491687153 rho 57: -0.09168717760514207 rho 60: -0.1239606605892366 rho 61: 0.13800141862943802 rho 69: 0.10269095820413282 rho 72: -0.09462958643919699 rho 82: -0.08002594037040578 rho 85: -0.08198325512480845 rho 93: 0.092203788733838 rho 109: 0.09582379018763891 rho 113: -0.09105599732819755 rho 119: 0.09171771966671013 rho 120: -0.13558724808031059
Se pensó que podía ser este modelo, pero es muy malo.
modelo=SARIMAX(zpre,
order=([1,5],0,[1,5]),
seasonal_order=([1],1,[1,2],12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
warn('Non-stationary starting autoregressive parameters'
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
warn('Non-invertible starting MA parameters found.'
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX([1, 5], 0, [1, 5])x(1, 1, 2, 12) | Log Likelihood | -5.099 |
| Date: | Sun, 27 Apr 2025 | AIC | 26.198 |
| Time: | 04:51:06 | BIC | 61.990 |
| Sample: | 0 | HQIC | 40.083 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 1.1476 | 32.468 | 0.035 | 0.972 | -62.489 | 64.784 |
| ar.L5 | -0.8775 | 38.486 | -0.023 | 0.982 | -76.308 | 74.553 |
| ma.L1 | -0.9317 | 23.880 | -0.039 | 0.969 | -47.736 | 45.872 |
| ma.L5 | 0.3314 | 36.742 | 0.009 | 0.993 | -71.681 | 72.344 |
| ar.S.L12 | -0.3697 | 133.649 | -0.003 | 0.998 | -262.317 | 261.577 |
| ma.S.L12 | -0.1630 | 41.007 | -0.004 | 0.997 | -80.536 | 80.210 |
| ma.S.L24 | -0.7558 | 42.337 | -0.018 | 0.986 | -83.735 | 82.223 |
| sigma2 | 0.4434 | 8.359 | 0.053 | 0.958 | -15.941 | 16.827 |
| Ljung-Box (L1) (Q): | 294.50 | Jarque-Bera (JB): | 2134861.83 |
|---|---|---|---|
| Prob(Q): | 0.00 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.00 | Skew: | 15.72 |
| Prob(H) (two-sided): | 0.00 | Kurtosis: | 282.43 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 6.39e+17. Standard errors may be unstable.
modelo.resid.mean()
-7.786026841214317e+48
También este, pero no es invertible
modelo=SARIMAX(zpre,
order=(0,0, [5]),
seasonal_order=(1,1,[1,2],12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting. self._init_dates(dates, freq) c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting. self._init_dates(dates, freq)
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(0, 0, [5])x(1, 1, [1, 2], 12) | Log Likelihood | -708.534 |
| Date: | Sun, 27 Apr 2025 | AIC | 1427.068 |
| Time: | 04:51:13 | BIC | 1449.437 |
| Sample: | 0 | HQIC | 1435.745 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ma.L5 | -0.1707 | 0.044 | -3.861 | 0.000 | -0.257 | -0.084 |
| ar.S.L12 | -0.8858 | 0.147 | -6.011 | 0.000 | -1.175 | -0.597 |
| ma.S.L12 | 0.2236 | 0.140 | 1.599 | 0.110 | -0.050 | 0.498 |
| ma.S.L24 | -0.6343 | 0.088 | -7.209 | 0.000 | -0.807 | -0.462 |
| sigma2 | 0.5151 | 0.019 | 27.358 | 0.000 | 0.478 | 0.552 |
| Ljung-Box (L1) (Q): | 10.39 | Jarque-Bera (JB): | 283.30 |
|---|---|---|---|
| Prob(Q): | 0.00 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | 0.66 |
| Prob(H) (two-sided): | 0.18 | Kurtosis: | 5.96 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
invertible(modelo)
El polinomio de la parte media móvil es: {'p5': -0.17066687867489047, 'p12': 0.2235561375559531, 'p24': -0.6342522946355272}
El grado del polinomio es: 24
Raíces del polinomio característico: [-1.03809773+0.j -0.97339645+0.26906311j -0.97339645-0.26906311j
-0.88525875+0.51576505j -0.88525875-0.51576505j -0.71259807+0.7035491j
-0.71259807-0.7035491j -0.52270681+0.89327539j -0.52270681-0.89327539j
-0.26081491+0.98110341j -0.26081491-0.98110341j 0.00733524+1.03148712j
0.00733524-1.03148712j 1.02333621+0.j 0.97341295+0.25216589j
0.97341295-0.25216589j 0.89999742+0.51572558j 0.89999742-0.51572558j
0.71256537+0.7204749j 0.71256537-0.7204749j 0.50801366+0.89331475j
0.50801366-0.89331475j 0.2608311 +0.96414954j 0.2608311 -0.96414954j]
Módulo de las raíces: [1.03809773 1.0098988 1.0098988 1.02454704 1.02454704 1.00138771
1.00138771 1.0349702 1.0349702 1.01517895 1.01517895 1.03151321
1.03151321 1.02333621 1.00554483 1.00554483 1.03728888 1.03728888
1.01332793 1.01332793 1.02766197 1.02766197 0.99880789 0.99880789]
¿Las raíces están fuera del círculo unitario? False
El modelo no es invertible
también este, pero no es invertible.
modelo=SARIMAX(zpre,
order=([12],0, [5]),
seasonal_order=(0,1,[1,2],12)).fit()
modelo.summary()
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting. self._init_dates(dates, freq) c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting. self._init_dates(dates, freq)
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX([12], 0, [5])x(0, 1, [1, 2], 12) | Log Likelihood | -708.534 |
| Date: | Sun, 27 Apr 2025 | AIC | 1427.068 |
| Time: | 04:51:22 | BIC | 1449.437 |
| Sample: | 0 | HQIC | 1435.745 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L12 | -0.8858 | 0.147 | -6.011 | 0.000 | -1.175 | -0.597 |
| ma.L5 | -0.1707 | 0.044 | -3.861 | 0.000 | -0.257 | -0.084 |
| ma.S.L12 | 0.2236 | 0.140 | 1.599 | 0.110 | -0.050 | 0.498 |
| ma.S.L24 | -0.6343 | 0.088 | -7.208 | 0.000 | -0.807 | -0.462 |
| sigma2 | 0.5151 | 0.019 | 27.358 | 0.000 | 0.478 | 0.552 |
| Ljung-Box (L1) (Q): | 10.39 | Jarque-Bera (JB): | 283.30 |
|---|---|---|---|
| Prob(Q): | 0.00 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | 0.66 |
| Prob(H) (two-sided): | 0.18 | Kurtosis: | 5.96 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
invertible(modelo)
El polinomio de la parte media móvil es: {'p5': -0.1706643002131486, 'p12': 0.22358740803256344, 'p24': -0.6342671085676491}
El grado del polinomio es: 24
Raíces del polinomio característico: [-1.03809804+0.j -0.97339404+0.26906235j -0.97339404-0.26906235j
-0.88525946+0.5157653j -0.88525946-0.5157653j -0.71259631+0.70354748j
-0.71259631-0.70354748j -0.52270685+0.89327583j -0.52270685-0.89327583j
-0.26081426+0.98110091j -0.26081426-0.98110091j 0.00733501+1.03148763j
0.00733501-1.03148763j 1.023337 +0.j 0.97341054+0.25216538j
0.97341054-0.25216538j 0.89999765+0.51572583j 0.89999765-0.51572583j
0.7125636 +0.72047304j 0.7125636 -0.72047304j 0.50801417+0.89331519j
0.50801417-0.89331519j 0.26083046+0.96414729j 0.26083046-0.96414729j]
Módulo de las raíces: [1.03809804 1.00989628 1.00989628 1.02454778 1.02454778 1.00138532
1.00138532 1.0349706 1.0349706 1.01517637 1.01517637 1.03151371
1.03151371 1.023337 1.00554237 1.00554237 1.03728921 1.03728921
1.01332536 1.01332536 1.0276626 1.0276626 0.99880555 0.99880555]
¿Las raíces están fuera del círculo unitario? False
El modelo no es invertible
MODELO PARA LA SERIE TRANSFORMADA¶
fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )
Valores de autocorrelacion significativos: r1: 0.1881206416574384 r2: 0.12601374062799692 r6: -0.10023577930193152 r12: -0.40245032162240263 r33: 0.09802046981427195 r45: -0.11985655705785987
facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])
Valores de autocorrelacion parcial significativos: rho 1: 0.18841139999075743 rho 2: 0.09425081653205467 rho 6: -0.0848567418000625 rho 12: -0.4319892618507997 rho 13: 0.08136260875048439 rho 17: -0.11385209887045368 rho 22: -0.08401445717662205 rho 24: -0.2809725408252081 rho 25: 0.10755279636113212 rho 30: -0.09673228306484968 rho 33: 0.1105360450418002 rho 36: -0.14602967562386657 rho 42: -0.11769525648448687 rho 46: -0.09378255882994088 rho 48: -0.16710968828081296 rho 50: 0.09151867173223789 rho 55: 0.09278816335947565 rho 60: -0.09971581325023068 rho 72: -0.13002571917786468 rho 78: -0.07964448430687596 rho 84: -0.07986957194068202 rho 92: -0.13966165907823416 rho 93: 0.11600081821072342 rho 94: -0.1221857229746919 rho 113: -0.08611614522897519 rho 118: -0.08199830520693245 rho 120: -0.11627799482364788
1 MODELO¶
Por la FAC y FACP se cree que puede ser el siguiente modelo:
modelo=SARIMAX(ypre,
order=([1, 2, 6],0,1),
seasonal_order=(1,1,6,12)).fit()
modelo.summary()
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX([1, 2, 6], 0, 1)x(1, 1, [1, 2, 3, 4, 5, 6], 12) | Log Likelihood | -617.741 |
| Date: | Sun, 27 Apr 2025 | AIC | 1259.482 |
| Time: | 18:48:15 | BIC | 1313.169 |
| Sample: | 0 | HQIC | 1280.309 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.3566 | 0.144 | 2.482 | 0.013 | 0.075 | 0.638 |
| ar.L2 | 0.0211 | 0.062 | 0.342 | 0.732 | -0.100 | 0.142 |
| ar.L6 | -0.2754 | 0.036 | -7.653 | 0.000 | -0.346 | -0.205 |
| ma.L1 | -0.1092 | 0.149 | -0.731 | 0.465 | -0.402 | 0.184 |
| ar.S.L12 | -0.8536 | 0.246 | -3.468 | 0.001 | -1.336 | -0.371 |
| ma.S.L12 | 0.1009 | 0.245 | 0.412 | 0.680 | -0.379 | 0.581 |
| ma.S.L24 | -0.7034 | 0.191 | -3.691 | 0.000 | -1.077 | -0.330 |
| ma.S.L36 | 0.0214 | 0.049 | 0.435 | 0.664 | -0.075 | 0.118 |
| ma.S.L48 | 0.0154 | 0.054 | 0.286 | 0.775 | -0.090 | 0.121 |
| ma.S.L60 | 0.0169 | 0.053 | 0.319 | 0.750 | -0.087 | 0.121 |
| ma.S.L72 | 0.0063 | 0.055 | 0.114 | 0.909 | -0.101 | 0.114 |
| sigma2 | 0.3874 | 0.020 | 19.120 | 0.000 | 0.348 | 0.427 |
| Ljung-Box (L1) (Q): | 0.02 | Jarque-Bera (JB): | 10.86 |
|---|---|---|---|
| Prob(Q): | 0.88 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | -0.20 |
| Prob(H) (two-sided): | 0.17 | Kurtosis: | 3.50 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
RESIDUOS INDEPENDIENTES¶
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 0.066845 | NaN |
| 2 | 1.537653 | NaN |
| 3 | 1.886894 | NaN |
| 4 | 4.744317 | NaN |
| 5 | 19.994274 | NaN |
| 6 | 20.002045 | NaN |
| 7 | 20.551066 | NaN |
| 8 | 20.563646 | NaN |
| 9 | 25.131429 | NaN |
| 10 | 25.236634 | NaN |
| 11 | 31.550176 | NaN |
| 12 | 31.652621 | 1.843656e-08 |
| 13 | 31.779184 | 1.256718e-07 |
| 14 | 31.817201 | 5.718884e-07 |
| 15 | 31.850384 | 2.052629e-06 |
| 16 | 34.932704 | 1.551903e-06 |
| 17 | 42.942269 | 1.197523e-07 |
| 18 | 48.759895 | 2.528439e-08 |
| 19 | 50.340543 | 3.514875e-08 |
| 20 | 50.340543 | 9.294961e-08 |
| 21 | 51.913020 | 1.184210e-07 |
| 22 | 52.744140 | 1.998096e-07 |
| 23 | 57.847442 | 5.554601e-08 |
| 24 | 57.847451 | 1.270340e-07 |
| 25 | 59.668596 | 1.341066e-07 |
| 26 | 59.677107 | 2.865712e-07 |
| 27 | 59.912799 | 5.413768e-07 |
| 28 | 60.290463 | 9.410050e-07 |
| 29 | 60.475407 | 1.713942e-06 |
| 30 | 69.299462 | 1.202209e-07 |
| 31 | 69.362452 | 2.313866e-07 |
| 32 | 69.561924 | 4.130052e-07 |
| 33 | 73.056696 | 2.166748e-07 |
| 34 | 74.446295 | 2.456685e-07 |
| 35 | 76.961663 | 1.840763e-07 |
| 36 | 76.961687 | 3.388905e-07 |
| 37 | 78.817354 | 3.185844e-07 |
| 38 | 80.532525 | 3.135244e-07 |
| 39 | 83.152029 | 2.243965e-07 |
| 40 | 83.267861 | 3.812489e-07 |
| 41 | 84.188435 | 4.853694e-07 |
| 42 | 86.097859 | 4.396425e-07 |
| 43 | 86.761902 | 6.023649e-07 |
| 44 | 86.884845 | 9.763344e-07 |
| 45 | 90.327760 | 5.298933e-07 |
| 46 | 91.250891 | 6.557028e-07 |
| 47 | 98.721616 | 9.467501e-08 |
| 48 | 98.935868 | 1.493370e-07 |
Hay dependencia en los residuos.
RESIDUOS CON MEDIA CERO¶
modelo.resid.mean()
-0.017262744929477162
ttest_1samp(modelo.resid, 0)
TtestResult(statistic=-0.6945529062973945, pvalue=0.4875803217531771, df=659)
La media de los residuos puede ser 0.
GRÁFICO DE RESIDUOS¶
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)
axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
line.set_linewidth(0.5)
# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas
# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
# line.set_markerfacecolor('green')
# line.set_markeredgecolor('white')
line.set_alpha(0.2)
# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)
plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
2 MODELO¶
Pero no captura completamente la dependencia temporal, por lo que se agregó la parte completa para la parte AR
modelo=SARIMAX(ypre,
order=(6,0,1),
seasonal_order=(1,1,1,12)).fit()
modelo.summary()
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(6, 0, 1)x(1, 1, 1, 12) | Log Likelihood | -602.344 |
| Date: | Sun, 27 Apr 2025 | AIC | 1224.688 |
| Time: | 18:50:32 | BIC | 1269.426 |
| Sample: | 0 | HQIC | 1242.043 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.4302 | 0.148 | 2.908 | 0.004 | 0.140 | 0.720 |
| ar.L2 | 0.0474 | 0.067 | 0.702 | 0.483 | -0.085 | 0.180 |
| ar.L3 | -0.0470 | 0.047 | -0.995 | 0.320 | -0.140 | 0.046 |
| ar.L4 | -0.0870 | 0.043 | -2.013 | 0.044 | -0.172 | -0.002 |
| ar.L5 | -0.1800 | 0.045 | -3.968 | 0.000 | -0.269 | -0.091 |
| ar.L6 | -0.2197 | 0.056 | -3.904 | 0.000 | -0.330 | -0.109 |
| ma.L1 | -0.1680 | 0.154 | -1.089 | 0.276 | -0.470 | 0.134 |
| ar.S.L12 | 0.0208 | 0.046 | 0.451 | 0.652 | -0.070 | 0.111 |
| ma.S.L12 | -0.8907 | 0.029 | -30.470 | 0.000 | -0.948 | -0.833 |
| sigma2 | 0.3657 | 0.019 | 19.149 | 0.000 | 0.328 | 0.403 |
| Ljung-Box (L1) (Q): | 0.02 | Jarque-Bera (JB): | 5.80 |
|---|---|---|---|
| Prob(Q): | 0.89 | Prob(JB): | 0.05 |
| Heteroskedasticity (H): | 0.81 | Skew: | -0.17 |
| Prob(H) (two-sided): | 0.12 | Kurtosis: | 3.32 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
RESIDUOS INDEPENDIENTES¶
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 0.078029 | NaN |
| 2 | 0.081745 | NaN |
| 3 | 0.110251 | NaN |
| 4 | 0.667736 | NaN |
| 5 | 1.809556 | NaN |
| 6 | 3.209615 | NaN |
| 7 | 3.289842 | NaN |
| 8 | 3.816187 | NaN |
| 9 | 6.473299 | NaN |
| 10 | 7.880155 | 0.004998 |
| 11 | 9.034590 | 0.010919 |
| 12 | 9.198592 | 0.026764 |
| 13 | 9.646686 | 0.046819 |
| 14 | 9.995740 | 0.075356 |
| 15 | 10.036647 | 0.123117 |
| 16 | 10.813321 | 0.146972 |
| 17 | 20.147979 | 0.009790 |
| 18 | 24.051508 | 0.004221 |
| 19 | 24.928482 | 0.005483 |
| 20 | 25.176289 | 0.008590 |
| 21 | 25.707594 | 0.011804 |
| 22 | 30.300176 | 0.004262 |
| 23 | 32.585774 | 0.003304 |
| 24 | 32.718036 | 0.005134 |
| 25 | 33.218651 | 0.006912 |
| 26 | 33.280245 | 0.010384 |
| 27 | 33.486398 | 0.014565 |
| 28 | 33.781554 | 0.019498 |
| 29 | 33.785406 | 0.027614 |
| 30 | 39.683798 | 0.008123 |
| 31 | 39.694743 | 0.011734 |
| 32 | 39.795172 | 0.016202 |
| 33 | 43.588079 | 0.008516 |
| 34 | 43.726951 | 0.011633 |
| 35 | 45.030235 | 0.011682 |
| 36 | 45.950944 | 0.012878 |
| 37 | 49.181595 | 0.007969 |
| 38 | 50.128815 | 0.008747 |
| 39 | 52.006371 | 0.007605 |
| 40 | 52.303949 | 0.009733 |
| 41 | 53.661655 | 0.009590 |
| 42 | 55.762429 | 0.007912 |
| 43 | 56.156343 | 0.009780 |
| 44 | 56.305459 | 0.012669 |
| 45 | 60.129401 | 0.007053 |
| 46 | 63.211636 | 0.004624 |
| 47 | 69.223074 | 0.001459 |
| 48 | 69.868377 | 0.001735 |
Y el modelo captura correctamente la dependencia temporal.
import modelo_admisible as modadm
modadm.estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.4302180830507938, 'p2': -0.04735397876460526, 'p3': 0.04698747454292404, 'p4': 0.08700716618181399, 'p5': 0.1799715794603173, 'p6': 0.21973455421960683, 'p12': -0.02082020047854175}
Raíces del polinomio característico: [-1.51585027+0.j -1.16433253+0.66675602j -1.16433253-0.66675602j
-0.70531602+1.36923083j -0.70531602-1.36923083j 0.00589534+1.27354789j
0.00589534-1.27354789j 0.89804478+1.3208691j 0.89804478-1.3208691j
1.62535851+0.j 0.91095431+0.5227982j 0.91095431-0.5227982j ]
Módulo de las raíces: [1.51585027 1.34172793 1.34172793 1.54021549 1.54021549 1.27356154
1.27356154 1.59724125 1.59724125 1.62535851 1.0503122 1.0503122 ]
¿Las raíces están fuera del círculo unitario? True
El modelo es estacionario
:)
modadm.invertible(modelo)
El polinomio de la parte media móvil es: {'p1': -0.16798341167047423, 'p12': -0.8906607240754797}
Raíces del polinomio característico: [-1.02312938+0.j -0.88152548+0.51636208j -0.88152548-0.51636208j
-0.49857917+0.88672189j -0.49857917-0.88672189j 0.01419601+1.01059688j
0.01419601-1.01059688j 0.51292531+0.86213479j 0.51292531-0.86213479j
0.99443323+0.j 0.86733141+0.49151162j 0.86733141-0.49151162j]
Módulo de las raíces: [1.02312938 1.02162467 1.02162467 1.01727917 1.01727917 1.01069658
1.01069658 1.00317933 1.00317933 0.99443323 0.99691898 0.99691898]
¿Las raíces están fuera del círculo unitario? False
El modelo no es invertible
3 MODELO¶
modelo=SARIMAX(ypre,
order=(6,0,0),
seasonal_order=(1,1,1,12)).fit()
modelo.summary()
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(6, 0, 0)x(1, 1, [1], 12) | Log Likelihood | -603.289 |
| Date: | Sun, 27 Apr 2025 | AIC | 1224.578 |
| Time: | 19:07:19 | BIC | 1264.843 |
| Sample: | 0 | HQIC | 1240.198 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.2721 | 0.034 | 7.929 | 0.000 | 0.205 | 0.339 |
| ar.L2 | 0.1022 | 0.041 | 2.473 | 0.013 | 0.021 | 0.183 |
| ar.L3 | -0.0255 | 0.040 | -0.644 | 0.520 | -0.103 | 0.052 |
| ar.L4 | -0.0900 | 0.041 | -2.215 | 0.027 | -0.170 | -0.010 |
| ar.L5 | -0.1974 | 0.039 | -5.011 | 0.000 | -0.275 | -0.120 |
| ar.L6 | -0.2604 | 0.038 | -6.796 | 0.000 | -0.336 | -0.185 |
| ar.S.L12 | 0.0196 | 0.046 | 0.426 | 0.670 | -0.071 | 0.110 |
| ma.S.L12 | -0.8851 | 0.029 | -30.327 | 0.000 | -0.942 | -0.828 |
| sigma2 | 0.3670 | 0.019 | 19.130 | 0.000 | 0.329 | 0.405 |
| Ljung-Box (L1) (Q): | 0.22 | Jarque-Bera (JB): | 6.43 |
|---|---|---|---|
| Prob(Q): | 0.64 | Prob(JB): | 0.04 |
| Heteroskedasticity (H): | 0.80 | Skew: | -0.19 |
| Prob(H) (two-sided): | 0.10 | Kurtosis: | 3.30 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
PRINCIPIO DE PARSIMONIA¶
parsimonia(modelo)
Hay coeficientes no significativos, no se cumple el principio de parsimonia
MODELO ADMISIBLE¶
modadm.estacionario(modelo)
El polinomio de la parte autoregresiva es: {'p1': -0.2721484850682655, 'p2': -0.10222956018637822, 'p3': 0.02552839449395747, 'p4': 0.09002860439631356, 'p5': 0.1974326132169234, 'p6': 0.2604147375399908, 'p12': -0.01958784345380709}
Raíces del polinomio característico: [-1.5324563 +0.j -0.72632726+1.38599378j -0.72632726-1.38599378j
-1.13728973+0.63130766j -1.13728973-0.63130766j 0.91165759+1.36955392j
0.91165759-1.36955392j -0.03290209+1.25835037j -0.03290209-1.25835037j
1.6720367 +0.j 0.91507128+0.53286918j 0.91507128-0.53286918j]
Módulo de las raíces: [1.5324563 1.56477795 1.56477795 1.30076027 1.30076027 1.64523479
1.64523479 1.25878045 1.25878045 1.6720367 1.0589169 1.0589169 ]
¿Las raíces están fuera del círculo unitario? True
El modelo es estacionario
:)
modadm.invertible(modelo)
El polinomio de la parte media móvil es: {'p12': -0.8851100185522583}
Raíces del polinomio característico: [-1.01022217e+00+0.j -8.74878063e-01+0.50511109j
-8.74878063e-01-0.50511109j -5.05111085e-01+0.87487806j
-5.05111085e-01-0.87487806j -5.55111512e-17+1.01022217j
-5.55111512e-17-1.01022217j 5.05111085e-01+0.87487806j
5.05111085e-01-0.87487806j 1.01022217e+00+0.j
8.74878063e-01+0.50511109j 8.74878063e-01-0.50511109j]
Módulo de las raíces: [1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217
1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217]
¿Las raíces están fuera del círculo unitario? True
El modelo es invertible
:)
RESIDUOS INDEPENDIENTES¶
acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)
| lb_stat | lb_pvalue | |
|---|---|---|
| 1 | 0.001254 | NaN |
| 2 | 0.127263 | NaN |
| 3 | 0.145069 | NaN |
| 4 | 0.355249 | NaN |
| 5 | 1.452317 | NaN |
| 6 | 3.552256 | NaN |
| 7 | 3.845680 | NaN |
| 8 | 4.082125 | NaN |
| 9 | 6.692169 | 0.009684 |
| 10 | 7.832130 | 0.019919 |
| 11 | 9.102569 | 0.027958 |
| 12 | 9.256310 | 0.055003 |
| 13 | 9.544643 | 0.089215 |
| 14 | 9.734154 | 0.136303 |
| 15 | 9.887684 | 0.195028 |
| 16 | 10.818577 | 0.212192 |
| 17 | 19.823382 | 0.019034 |
| 18 | 23.899218 | 0.007872 |
| 19 | 24.703424 | 0.010073 |
| 20 | 25.013777 | 0.014758 |
| 21 | 25.631668 | 0.019043 |
| 22 | 29.998524 | 0.007635 |
| 23 | 32.696622 | 0.005169 |
| 24 | 32.798350 | 0.007858 |
| 25 | 33.458538 | 0.009854 |
| 26 | 33.506221 | 0.014484 |
| 27 | 33.689446 | 0.019989 |
| 28 | 33.907394 | 0.026758 |
| 29 | 33.908062 | 0.037069 |
| 30 | 40.197221 | 0.010252 |
| 31 | 40.244748 | 0.014426 |
| 32 | 40.383594 | 0.019444 |
| 33 | 44.092471 | 0.010590 |
| 34 | 44.309675 | 0.013997 |
| 35 | 45.634276 | 0.013925 |
| 36 | 46.439778 | 0.015687 |
| 37 | 49.474952 | 0.010282 |
| 38 | 50.724475 | 0.010416 |
| 39 | 52.415648 | 0.009474 |
| 40 | 52.583380 | 0.012372 |
| 41 | 53.835941 | 0.012450 |
| 42 | 55.965773 | 0.010224 |
| 43 | 56.386506 | 0.012439 |
| 44 | 56.652499 | 0.015534 |
| 45 | 60.570646 | 0.008572 |
| 46 | 63.146117 | 0.006367 |
| 47 | 69.344542 | 0.001975 |
| 48 | 69.962021 | 0.002346 |
Hay dependencia en los primeros dos residuos.
RESIDUOS CON MEDIA CERO¶
modelo.resid.mean()
-0.03293748637459734
ttest_1samp(modelo.resid, 0)
TtestResult(statistic=-1.3520655710406162, pvalue=0.1768181651558442, df=659)
La media de los residuos puede ser 0.
RESIDUOS CON VARIANZA CONSTANTE¶
het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))
(8.59265039527976, 0.0033752241057990228, 8.679613399395905, 0.003331636731645561)
Los residuos tienen varianza constante.
RESIDUOS CON DISTRIBUCIÓN NORMAL¶
print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))
SignificanceResult(statistic=6.693009944604184, pvalue=0.03520718942157648) (0.02689680139562356, 0.374813270764381)
Los residuos puede que sigan una distribución normal.
normal_relajacion(modelo.resid)
70.00% de los residuos están dentro de ±1σ (esperado ≈ 68%) 95.30% de los residuos están dentro de ±2σ (esperado ≈ 95%) 99.55% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)
Pero son bastante parecidos.
GRÁFICO DE RESIDUOS¶
fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)
axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
line.set_linewidth(0.5)
# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas
# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
# line.set_markerfacecolor('green')
# line.set_markeredgecolor('white')
line.set_alpha(0.2)
# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)
plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()
MODELO SELECCIONADO¶
Para la serie sin transformar, este es el modelo seleccionado:
$$\text{ARIMA}(13, 0, 0)(0, 1, 4)_{12}$$Nos quedamos con este modelo porque es el que tiene los p-valores más altos en la prueba de independencia de los residuos. Cumple todos los supuestos, menos el principio de parsimonia, y sus residuos siguen una distribución normal relajada. Tiene un $\text{AIC} = 1367.369$
Para la serie transformada, este es el modelo seleccionado:
$$\text{ARIMA}(6, 0, 0)(1, 1, 1)_{12}$$Cumple todos los supuestos, menos el principio de parsimonia, y tiene un $\text{AIC} = 1224.578$
modelo = SARIMAX(zpre,
order=(13, 0, 0),
seasonal_order=(0, 1, 4, 12)).fit()
modelo.summary()
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) | Log Likelihood | -665.247 |
| Date: | Sun, 27 Apr 2025 | AIC | 1366.493 |
| Time: | 19:21:59 | BIC | 1447.023 |
| Sample: | 0 | HQIC | 1397.733 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.1420 | 0.035 | 4.086 | 0.000 | 0.074 | 0.210 |
| ar.L2 | -0.0189 | 0.037 | -0.516 | 0.606 | -0.090 | 0.053 |
| ar.L3 | 0.0496 | 0.037 | 1.349 | 0.177 | -0.022 | 0.122 |
| ar.L4 | -0.0694 | 0.042 | -1.646 | 0.100 | -0.152 | 0.013 |
| ar.L5 | -0.1131 | 0.050 | -2.274 | 0.023 | -0.211 | -0.016 |
| ar.L6 | -0.1124 | 0.060 | -1.872 | 0.061 | -0.230 | 0.005 |
| ar.L7 | -0.0152 | 0.051 | -0.299 | 0.765 | -0.115 | 0.085 |
| ar.L8 | -0.0505 | 0.044 | -1.154 | 0.249 | -0.136 | 0.035 |
| ar.L9 | 0.0363 | 0.035 | 1.030 | 0.303 | -0.033 | 0.105 |
| ar.L10 | -0.0074 | 0.033 | -0.225 | 0.822 | -0.072 | 0.057 |
| ar.L11 | 0.0991 | 0.037 | 2.677 | 0.007 | 0.027 | 0.172 |
| ar.L12 | 0.5353 | 0.082 | 6.554 | 0.000 | 0.375 | 0.695 |
| ar.L13 | -0.0519 | 0.042 | -1.247 | 0.213 | -0.133 | 0.030 |
| ma.S.L12 | -1.3733 | 0.097 | -14.094 | 0.000 | -1.564 | -1.182 |
| ma.S.L24 | 0.3639 | 0.086 | 4.217 | 0.000 | 0.195 | 0.533 |
| ma.S.L36 | 0.0957 | 0.060 | 1.599 | 0.110 | -0.022 | 0.213 |
| ma.S.L48 | -0.0740 | 0.039 | -1.888 | 0.059 | -0.151 | 0.003 |
| sigma2 | 0.4357 | 0.028 | 15.586 | 0.000 | 0.381 | 0.491 |
| Ljung-Box (L1) (Q): | 0.11 | Jarque-Bera (JB): | 344.46 |
|---|---|---|---|
| Prob(Q): | 0.74 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | 0.87 |
| Prob(H) (two-sided): | 0.18 | Kurtosis: | 6.12 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
modelo2 = SARIMAX(ypre,
order=(6, 0, 0),
seasonal_order=(1, 1, 1, 12)).fit()
modelo2.summary()
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(6, 0, 0)x(1, 1, [1], 12) | Log Likelihood | -603.289 |
| Date: | Sun, 27 Apr 2025 | AIC | 1224.578 |
| Time: | 19:20:28 | BIC | 1264.843 |
| Sample: | 0 | HQIC | 1240.198 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.2721 | 0.034 | 7.929 | 0.000 | 0.205 | 0.339 |
| ar.L2 | 0.1022 | 0.041 | 2.473 | 0.013 | 0.021 | 0.183 |
| ar.L3 | -0.0255 | 0.040 | -0.644 | 0.520 | -0.103 | 0.052 |
| ar.L4 | -0.0900 | 0.041 | -2.215 | 0.027 | -0.170 | -0.010 |
| ar.L5 | -0.1974 | 0.039 | -5.011 | 0.000 | -0.275 | -0.120 |
| ar.L6 | -0.2604 | 0.038 | -6.796 | 0.000 | -0.336 | -0.185 |
| ar.S.L12 | 0.0196 | 0.046 | 0.426 | 0.670 | -0.071 | 0.110 |
| ma.S.L12 | -0.8851 | 0.029 | -30.327 | 0.000 | -0.942 | -0.828 |
| sigma2 | 0.3670 | 0.019 | 19.130 | 0.000 | 0.329 | 0.405 |
| Ljung-Box (L1) (Q): | 0.22 | Jarque-Bera (JB): | 6.43 |
|---|---|---|---|
| Prob(Q): | 0.64 | Prob(JB): | 0.04 |
| Heteroskedasticity (H): | 0.80 | Skew: | -0.19 |
| Prob(H) (two-sided): | 0.10 | Kurtosis: | 3.30 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
MÉTODO FORECAST¶
# Pronóstico para 12 meses hacia el futuro
forecast = modelo.forecast(steps=12)
forecast2 = modelo2.forecast(steps=12)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`. return get_prediction_index( c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception. return get_prediction_index(
# Obtener pronóstico con intervalos
pred = modelo.get_forecast(steps=12)
pred2 = modelo2.get_forecast(steps=12)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`. return get_prediction_index( c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception. return get_prediction_index(
# Media pronosticada
forecast_mean = pred.predicted_mean
forecast_mean
660 -0.276578 661 -0.661192 662 -0.642475 663 -0.132995 664 0.046452 665 0.706953 666 0.536280 667 0.543615 668 0.539949 669 -0.053797 670 -0.541525 671 -0.504177 Name: predicted_mean, dtype: float64
# Media pronosticada
forecast2_mean = pred2.predicted_mean
forecast2_mean # El output es un array de numpy, porque tiene la transformacion
array([-0.8142261 , -0.58943894, -0.67426909, -0.42378983, 0.14539092,
0.48758714, 0.50076814, 0.57193519, 0.59488235, 0.21526718,
-0.35086441, -0.49899387])
# aplicar la inversa de la estandarización
# Calcular media y desviación estándar de la serie original
mu = np.mean(pre)
sigma = np.std(pre, ddof=0) # ddof=0 porque así lo usa scipy.stats.zscore por default
# Desestandarizar
forecast_mean = forecast_mean * sigma + mu
forecast_mean
660 540.502158 661 201.870616 662 218.349640 663 666.919279 664 824.911936 665 1406.446553 666 1256.178392 667 1262.636919 668 1259.408611 669 736.647934 670 307.230864 671 340.113522 Name: predicted_mean, dtype: float64
forecast2_mean = forecast2_mean * sigma + mu
forecast2_mean
array([ 67.13264749, 265.04529754, 190.35704854, 410.89015873,
912.02229063, 1213.30710951, 1224.91224955, 1287.57089294,
1307.77459904, 973.54446573, 475.09695852, 344.67717081])
# Intervalos de confianza
int_conf = pred.conf_int()
int_conf
| lower y | upper y | |
|---|---|---|
| 660 | -1.566687 | 1.023425 |
| 661 | -1.969017 | 0.645214 |
| 662 | -1.962076 | 0.652155 |
| 663 | -1.451584 | 1.165524 |
| 664 | -1.273046 | 1.348105 |
| 665 | -0.629675 | 2.012863 |
| 666 | -0.807778 | 1.859318 |
| 667 | -0.787023 | 1.882975 |
| 668 | -0.791284 | 1.883888 |
| 669 | -1.390333 | 1.285407 |
| 670 | -1.872244 | 0.804378 |
| 671 | -1.848288 | 0.850488 |
# Desestandarizar los intervalos de confianza
int_conf = int_conf * sigma + mu
int_conf
| lower y | upper y | |
|---|---|---|
| 660 | -595.367749 | 1685.082779 |
| 661 | -949.596462 | 1352.088860 |
| 662 | -943.485554 | 1358.199735 |
| 663 | -494.026019 | 1810.193098 |
| 664 | -336.832977 | 1970.945813 |
| 665 | 229.619785 | 2556.228270 |
| 666 | 72.809546 | 2421.040094 |
| 667 | 91.083533 | 2441.868542 |
| 668 | 87.331818 | 2442.672581 |
| 669 | -440.097428 | 1915.743419 |
| 670 | -864.393291 | 1492.223637 |
| 671 | -843.301917 | 1532.821294 |
# Asignar los índices de fecha a los intervalos de confianza
forecast_index = pd.date_range(start=pre.index[-2], periods=13, freq='ME')[1:]
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3)) # Pone cada 3 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) # Formato de fecha
# Serie original
plt.plot(pre_total, label='Observado')
# Pronóstico
plt.plot(forecast_index, forecast_mean, label='Pronóstico', color=sns.color_palette()[1])
# Intervalos
plt.fill_between(forecast_index,
int_conf.iloc[:, 0], int_conf.iloc[:, 1],
color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.legend()
plt.title('Método Forecast para modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
plt.savefig('imagenes/05-01-metodo-forecast.svg', bbox_inches='tight')
plt.show()
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4)) # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y')) # Formato de fecha
# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']
# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, forecast_mean , label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
# Intervalos
plt.fill_between(forecast_index,
int_conf.iloc[:, 0], int_conf.iloc[:, 1],
color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-02-metodo-forecast-zoom.svg', bbox_inches='tight')
plt.show()
PRONÓSTICO ÓPTIMO¶
Se calcula en 04-Calculo-del-Pronostico-Optimo
pre_total.to_csv('04-01-Serie-Original.csv', index=False) # Guardar la serie original
modelo.resid.to_csv('04-02-Serie-Residuos.csv', index=False) # Guardar la serie de residuos
modelo.summary()
| Dep. Variable: | y | No. Observations: | 660 |
|---|---|---|---|
| Model: | SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12) | Log Likelihood | -665.685 |
| Date: | Sun, 27 Apr 2025 | AIC | 1367.369 |
| Time: | 04:52:23 | BIC | 1447.899 |
| Sample: | 0 | HQIC | 1398.609 |
| - 660 | |||
| Covariance Type: | opg |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| ar.L1 | 0.1369 | 0.035 | 3.963 | 0.000 | 0.069 | 0.205 |
| ar.L2 | -0.0190 | 0.036 | -0.522 | 0.602 | -0.090 | 0.052 |
| ar.L3 | 0.0501 | 0.037 | 1.363 | 0.173 | -0.022 | 0.122 |
| ar.L4 | -0.0697 | 0.042 | -1.662 | 0.097 | -0.152 | 0.013 |
| ar.L5 | -0.1117 | 0.050 | -2.254 | 0.024 | -0.209 | -0.015 |
| ar.L6 | -0.1101 | 0.060 | -1.839 | 0.066 | -0.227 | 0.007 |
| ar.L7 | -0.0104 | 0.051 | -0.205 | 0.838 | -0.110 | 0.089 |
| ar.L8 | -0.0511 | 0.044 | -1.170 | 0.242 | -0.137 | 0.034 |
| ar.L9 | 0.0331 | 0.035 | 0.944 | 0.345 | -0.036 | 0.102 |
| ar.L10 | -0.0093 | 0.033 | -0.281 | 0.778 | -0.074 | 0.055 |
| ar.L11 | 0.1034 | 0.037 | 2.802 | 0.005 | 0.031 | 0.176 |
| ar.L12 | 0.5385 | 0.080 | 6.699 | 0.000 | 0.381 | 0.696 |
| ar.L13 | -0.0460 | 0.041 | -1.115 | 0.265 | -0.127 | 0.035 |
| ma.S.L12 | -1.3756 | 0.100 | -13.725 | 0.000 | -1.572 | -1.179 |
| ma.S.L24 | 0.3674 | 0.086 | 4.275 | 0.000 | 0.199 | 0.536 |
| ma.S.L36 | 0.0973 | 0.060 | 1.635 | 0.102 | -0.019 | 0.214 |
| ma.S.L48 | -0.0784 | 0.039 | -1.994 | 0.046 | -0.155 | -0.001 |
| sigma2 | 0.4340 | 0.030 | 14.503 | 0.000 | 0.375 | 0.493 |
| Ljung-Box (L1) (Q): | 0.04 | Jarque-Bera (JB): | 347.08 |
|---|---|---|---|
| Prob(Q): | 0.85 | Prob(JB): | 0.00 |
| Heteroskedasticity (H): | 0.83 | Skew: | 0.87 |
| Prob(H) (two-sided): | 0.17 | Kurtosis: | 6.14 |
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
modelo.params
ar.L1 0.136917 ar.L2 -0.019014 ar.L3 0.050067 ar.L4 -0.069698 ar.L5 -0.111714 ar.L6 -0.110100 ar.L7 -0.010426 ar.L8 -0.051137 ar.L9 0.033088 ar.L10 -0.009268 ar.L11 0.103359 ar.L12 0.538506 ar.L13 -0.046024 ma.S.L12 -1.375601 ma.S.L24 0.367441 ma.S.L36 0.097350 ma.S.L48 -0.078355 sigma2 0.434001 dtype: float64
con coeficientes:
- $\phi_1 = -0.136917$
- $\phi_2 = 0.019014$
- $\phi_3 = -0.050067$
- $\phi_4 = 0.069698$
- $\phi_5 = 0.111714$
- $\phi_6 = 0.110100$
- $\phi_7 = 0.010426$
- $\phi_8 = 0.051137$
- $\phi_9 = -0.033088$
- $\phi_{10} = 0.009268$
- $\phi_{11} = -0.103359$
- $\phi_{12} = -0.538506$
- $\phi_{13} = 0.046024$
y para la parte estacional MA:
- $\Theta_1 = 1.375601$
- $\Theta_2 = -0.367441$
- $\Theta_3 = -0.097350$
- $\Theta_4 = 0.078355$
Donde $\varepsilon_t \sim \mathcal{N}(0, 1)$
Desarrollo para el pronóstico óptimo¶
Del lado izquierdo:
$$ (1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_{13} B^{13})(X_t - X_{t-12}) = X_t - X_{t-12} - \phi_1 (X_{t-1} - X_{t-13}) - \phi_2 (X_{t-2} - X_{t-14}) - \cdots - \phi_{13}(X_{t-13} - X_{t-25}) $$$$ = X_t - X_{t-12} - \sum_{h=1}^{13} \phi_k (X_{t-h} - X_{t-h-12}) $$Del lado derecho:
$$ (1 - \Theta_1 B^{12} - \Theta_2 B^{24} - \Theta_3 B^{36} - \Theta_4 B^{48})\varepsilon_t = \varepsilon_t - \Theta_1 \varepsilon_{t-12} - \Theta_2 \varepsilon_{t-24} - \Theta_3 \varepsilon_{t-36} - \Theta_4 \varepsilon_{t-48} $$despejando $X_t$:
$$ X_t = X_{t-12} + \sum_{k=1}^{13} \phi_k (X_{t-k} - X_{t-k-12}) + \varepsilon_t - \Theta_1 \varepsilon_{t-12} - \Theta_2 \varepsilon_{t-24} - \Theta_3 \varepsilon_{t-36} - \Theta_4 \varepsilon_{t-48} $$Por lo que, el pronóstico óptimo sería:
$$ X_t(h) = \underset{t}{\mathbb{E}}[X_{t+h-12}] + \sum_{k=1}^{13} \phi_k \left( \underset{t}{\mathbb{E}}[X_{t+h-k}] - \underset{t}{\mathbb{E}}[X_{t+h-k-12}] \right) + \underset{t}{\mathbb{E}}[\varepsilon_{t+h}] - \sum_{j=1}^4 \Theta_j \, \underset{t}{\mathbb{E}}[\varepsilon_{t+h-12j}] $$# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)
# Calcula de las espis
psi0 = -1
# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i>q
# pronostico optimo 04-Calculo-del-Pronostico-Optimo
# 893.2021421
# 218.8353082
# 14.55589463
# 46.06671573
# 102.6854446
# 909.4688428
# 328.2252003
# 431.4633135
# 1179.309604
# 567.0569539
# 460.930121
# 514.9098934
# Guardarlos en una serie
pronostico_optimo = pd.Series([893.2021421, 218.8353082, 14.55589463, 46.06671573, 102.6854446, 909.4688428, 328.2252003, 431.4633135, 1179.309604, 567.0569539, 460.930121, 514.9098934],
index=pd.date_range(start='2009-01-01', periods=12, freq='ME'))
# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)
# Calcula de las espis
psi0 = -1
# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i>q
phi1 = modelo.params[0]
teta12 = modelo.params[1]
psi1 = phi1*psi0
psi2 = phi1*psi1
psi3 = phi1*psi2
psi4 = phi1*psi3
psi5 = phi1*psi4
psi6 = phi1*psi5
psi7 = phi1*psi6
psi8 = phi1*psi7
psi9 = phi1*psi8
psi10 = phi1*psi9
psi11 = phi1*psi10
psi12 = phi1*psi11 + teta12
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:10: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` phi1 = modelo.params[0] C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:11: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` teta12 = modelo.params[1]
psis = [psi0, psi1, psi2, psi3, psi4, psi5, psi6, psi7, psi8, psi9, psi10, psi11, psi12]
sigmax = np.std(modelo.resid, ddof=0) # ddof=0 porque así lo usa scipy.stats.zscore por default
# 1.96 es el valor crítico para un intervalo de confianza del 95%
# primer intervalo
u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# segundo intervalo
u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# tercer intervalo
u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# cuarto intervalo
u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# quinto intervalo
u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# sexto intervalo
u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# septimo intervalo
u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# octavo intervalo
u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# noveno intervalo
u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# decimo intervalo
u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# decimo primero intervalo
u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# decimo segundo intervalo
u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:4: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:5: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:8: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:9: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:12: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:13: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:16: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:17: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:20: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:21: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:24: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:25: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:28: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:29: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:32: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:33: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:36: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:37: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:40: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:41: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:44: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:45: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos) C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:49: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
# Guardar los resultados en un dataframe
pronostico_optimo_df = pd.DataFrame({
'Pronóstico': pronostico_optimo,
'Límite Inferior': [l0, l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11],
'Límite Superior': [u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11]
}, index=pronostico_optimo.index)
import matplotlib.dates as mdates
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4)) # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y')) # Formato de fecha
# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']
# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
# Intervalos
plt.fill_between(forecast_index,
pronostico_optimo_df['Límite Inferior'], pronostico_optimo_df['Límite Superior'],
color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo con Esperanzas Condicionales')
plt.ylim(-1000, 6000)
plt.grid(True)
plt.savefig('imagenes/05-03-pronostico-optimo.svg', bbox_inches='tight')
plt.show()
COMPARACIÓN¶
# Comparar el pronóstico óptimo con el pronóstico del modelo SARIMAX
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4)) # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y')) # Formato de fecha
# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']
# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico Óptimo', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean, label='Método Forecast', color=sns.color_palette()[2], linewidth=3, alpha=0.7)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo vs Método Forecast')
plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-04-comparacion.svg', bbox_inches='tight')
plt.show()
ACTUALIZACIÓN DE PRONÓSTICO¶
# Conseguimos la primera observación de la serie test
valornuevo = pre_test[0]
# Estandarizamos
znv = (valornuevo - mu) / sigma
# Crear una nueva marca de tiempo para la nueva observación (usando el siguiente período después de la última fecha de la serie estandarizada)
new_date = zpre.index[-1] + pd.DateOffset(months=1)
new_obs = pd.Series(znv, index=[new_date])
# Actualizar la serie estandarizada con la nueva observación (si es necesario) utilizando pd.concat
zpre_updated = pd.concat([zpre, new_obs])
# Actualizar el modelo con la nueva observación sin re-calibrarlo
modelo_updated = modelo.append([new_obs], refit=False)
# Generar pronóstico para los próximos 12 meses (en escala estandarizada)
forecast_std_updated = modelo_updated.forecast(steps=11)
# Transformar el pronóstico de vuelta a la escala original
forecast_updated = forecast_std_updated * sigma + mu
# Asignar los índices de fecha a los intervalos de confianza
forecast_index2 = pd.date_range(start=pre.index[-1], periods=12, freq='ME')[1:]
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\2112643297.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` valornuevo = pre_test[0] c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`. return get_prediction_index( c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception. return get_prediction_index( c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`. return get_prediction_index(
plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4)) # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y')) # Formato de fecha
# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']
# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index2, forecast_updated , label='Pronóstico Actualizado', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean , label='Pronóstico Original', color=sns.color_palette()[2], linewidth=3, alpha=0.7)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/06-01-Actualizacion-de-pronostico-metodo-forecast.svg', bbox_inches='tight')
plt.show()